Imagine a mathematician wants to create this function:
$$F(s)=\sum_{n=0}^{\infty} \frac{1}{(n^2+1)^s}$$
Where would he start?
More precisely, how would he know the exact values of $F(s)$ for $s> 0$? And how would he extend this result for every natural $s$? And for real $s$?
Update:
If one may recall the following series:
$$\sum_{n=0}^\infty\frac1{n^2+x}=\frac1{2x}+\frac\pi{2\sqrt x\tanh(\pi\sqrt x)}$$
Then differentiating both sides $k$ times yields
$$\sum_{n=0}^\infty\frac1{(n^2+x)^{k+1}}=\frac1{2x^{k+1}}+\frac{(-1)^k\pi}{2\times k!}\frac{d^k}{dx^k}\frac1{\sqrt x\tanh(\pi\sqrt x)}$$
and evaluation at $x=1$ gives closed forms.
Update:
Using general Leibniz rule, we have
$$\frac{d^k}{dx^k}\frac{\coth(\pi\sqrt x)}{\sqrt x}=\sum_{p=0}^k\binom kp\frac{\Gamma(3/2+n-p)}{-2\sqrt\pi}(-1)^{n-p}x^{-1/2-n+p}(\coth(\pi\sqrt x))^{(p)}$$
We can handle the chain rule with the $n$th derivative of $\coth$ using Faà di Bruno's formula,
$$\small(\coth(\pi\sqrt x))^{(p)}=\sum_{q=1}^n\coth^{(q)}(\pi\sqrt x)B_{p,q}\left(\pi\frac12x^{-1/2},-\pi\frac14x^{-1/2},\dots,\sqrt\pi\frac{\Gamma(3/2+q)}{-4}(-1)^qx^{-1/2-q}\right)$$
where $B_{n,k}$ is the bell polynomial. The $q$th derivative of $\coth$ is then given when $q\ge1$:
$$\coth^{(q)}(x)=2^q(\coth(x)-1)\sum_{r=0}^q\frac{(-1)^rr!S_q^{(r)}}{2^r}(\coth(x)+1)^r$$
Putting all of this together,
$$\tiny F(k+1)=\frac12+\frac{(-1)^k\pi}{2(k!)}\left(\binom kp\frac{\Gamma(3/2+k)}{-2\sqrt\pi}(-1)^k(\coth(\pi)+\sum_{p=1}^k\sum_{q=1}^n\sum_{r=0}^q\binom kp\frac{\Gamma(3/2+k-p)}{-2\sqrt\pi}(-1)^{k+p+r}2^{q-r}(\coth(\pi)-1)(-1)^rr!S_q^{(r)}(\coth(\pi)+1)^r(\pi)B_{p,q}\left(\pi,-\pi\frac14,\dots,\sqrt\pi\frac{\Gamma(3/2+q)}{-4}(-1)^q\right)\right)$$
Old:
$$\frac12[\vartheta_3(0,r)+1]=\sum_{n=0}^\infty r^{n^2}$$
Term by term integration then reveals that
$$\frac1x\int_0^x\frac12[\vartheta_3(0,r)+1]\ dr=\sum_{n=0}^\infty\frac{x^{n^2}}{1+n^2}$$
Repeat this process over and over to get
$$\frac1{x_1}\int_0^{x_1}dx_2\frac1{x_2}\int_0^{x_2}dx_3\dots dx_k\frac1{x_k}\int_0^{x_k}dr\frac12[\vartheta_3(0,r)+1]=\sum_{n=0}^\infty\frac{x_1^{n^2}}{(1+n^2)^k}$$
And of course, evaluate at $x_1=1$ for your function. I do not see a closed form coming out of this, but it may be useful for discerning certain things about your function.
I do note, however, that in the case of $k=1$, the solution is given in this post:
$$\sum_{n=0}^\infty\frac1{1+n^2}=\int_0^1\frac12[\vartheta_3(0,r)+1]\ dr=\frac\pi{2\tanh\pi}+\frac12$$
Accepting, for argument's sake, that it is possible to create a function (the opposing view is that one discovers mathematical objects, rather than creating them), you have created the function simply by writing down that formula. Well, you have created it for real part of $s$ exceeding one-half, which is where the series converges. Knowing "exact values" (again, a term we can argue over, but I'll take it to mean a finite expression in terms of well-known constants, and leave it at that), well, that's generally very difficult. No one knows an exact value for $\zeta(3)$, for example, and for your function I suspect no one knows an exact value for $F(2)$. And as to extending your function to be defined for all complex $s$ (outside of a pole at $s=1/2$), that's what analytic continuation is for, so now you have a keyphrase to search for.
It doesn't have an Euler product.
Expanding $(1+n^{-2})^{-s} =\sum_{k=0}^\infty {-s \choose k} n^{-2k}$ you get $$F(s) =\sum_{n=0}^\infty (n^2+1)^{-s} = 1+2^{-s}+\sum_{k=0}^\infty {-s \choose k} (\zeta(2s+2k)-1)$$ so it is meromorphic on the whole complex plane, with poles at $s = \frac{1}{2}-k, k \in \mathbb{N}$
Using $\Gamma(s) a^{-s} =\int_0^\infty x^{s-1}e^{-ax}dx$ and $\theta(x) = \sum_{n=1}^\infty e^{-x n^2}$ you have $$G(s) = (F(s)-1-2^{-s}) \Gamma(s) = \int_0^\infty x^{s-1}e^{-x} \theta(x)dx, \qquad \Gamma(s) \zeta(2s) = \int_0^\infty x^{s-1}\theta(x) dx $$ where $\Gamma(s) \zeta(2s)-\frac{\sqrt{\pi}}{2(s-1/2)}+\frac{1}{2s}$ is entire.
From these poles location, we can deduce that for arbitrary large $N$ :
$\theta(x) = \frac{\sqrt{\pi}}{2}x^{-1/2}-\frac{1}{2}+o(x^N)$ as $x \to 0$ and hence $e^{-x}\theta(x) = \sum_{k\ge 0} \frac{x^k}{k!}(\frac{\sqrt{\pi}}{2}x^{-1/2}-\frac{1}{2})+o(x^N)$ as $x \to 0$
so that $$\lim_{s \to 1/2-k}(s+1/2-k)G(s)= \frac{\sqrt{\pi}}{2k!}, \qquad \lim_{s \to -k}(s-k)G(s)= \frac{-1/2}{(k+1)!}$$
Perhaps look into the theory of Dirichlet series? Your function is the series corresponding to the sequence $n \mapsto [\exists x\in \mathbb{N}(n=x^2+1)],$ where $[-]$ is the Iverson Bracket (google it.) The first few terms of the sequence are:
$$0,1,1,0,0,1\ldots$$
