# Constructing a matrix of solutions using findroot

I would like to solve the following equation numerically:

 Sinh[y]-(y+Tanh[x])Cosh[y]-y Sech[x]^2==0


I would like to do it by creating a list of x values:

xList = Array[# &, 10000, {-5, 5}];


and then find the root of v(x) for every x and store it in a list as well.

FindRoot[Sinh[y]-(y+Tanh[x])Cosh[y]-y Sech[x]^2==0,{v,vp}];


any ideas how to do it?

If you want to solve your equation only, just do it with FindInstance:

f = Sinh[y] - (y + Tanh[x]) Cosh[y] - y Sech[x]^2;
FindInstance[f == 0, {x, y}, Reals]

{{x -> 0, y -> 0}}


If you want to plot the Contour-Line F(x ,y(x)):

xContour = Range[-5, 5, 0.1];
yContour =
y /. FindInstance[(Sinh[y] - (y + Tanh[#]) Cosh[y] - y Sech[#]^2) == 0, {y}, Reals] & /@ xContour // N // Flatten;


This can be done also with Mathematicas ContourPlot.

ContourPlot[f == 0, {x, -5, 5}, {y, -3, 3}, GridLines -> Automatic]


Each point of the contour is a solution of the equation.

We wish to find $y(x)$ such that $f(x,y(x)) = 0$. We know that $f(0,0) = 0$; thus, $y(0) = 0$. We also know that if $f(x,y(x)) = 0$, then $$\frac{d}{dx} f(x,y(x)) = 0$$ as well. Thus, we can replace the condition $f(x,y(x)) = 0$ with the conditions $$\frac{d}{dx} f(x,y(x)) = 0, \qquad y(0) = 0.$$ These two equations are a first-order differential equation in $y(x)$, and an "initial condition" for $y(x)$. It's still not exactly soluble, but NDSolve can handle it numerically and return an InterpolatingFunction solution. You can then feed your desired values of x into the InterpolatingFunction using Map to get a corresponding table of $y$-values.

Implemented:

f[x_, y_] = Sinh[y] - (y + Tanh[x]) Cosh[y] - y Sech[x]^2;
soln = NDSolve[{y[0] == 0, D[f[x, y[x]], x] == 0}, y, {x, -5, 5}];
xList = Array[# &, 10000, {-5, 5}];
yList = First[y /. soln] /@ xList


This method has the advantage of not calling a numerical root-finder such as FindRoot or Solve 10,000 times, and so is pretty fast. On my machine, I get a result within about 40 milliseconds.

If you want to plot it, you can do that too:

Plot[y[x] /. soln, {x, -5, 5}]


Note that this method relies on three facts:

1. there was only one connected contour of $f(x,y) = 0$,
2. it led to a well-defined function $y(x)$, as opposed to multiple $y$-values corresponding to a single $x$, and
3. we knew a point that it passed through, namely $(0,0)$.

If you picked another $f$ you might not have such luck. Condition 3 would be easy enough to relax (a FindInstance call would probably give you a starting point for your initial conditions), but Conditions 1 and 2 would be harder to get around.

soln = ({N[#], y} /. FindRoot[
Sinh[y] - (y + Tanh[#]) Cosh[y] - y Sech[#]^2 == 0,
{y, -Sign[#]},
WorkingPrecision -> 20]) & /@
Range[-5, 5, 1/1000];


Looking at several examples

soln[[4995 ;; 5005]]

(*  {{-0.006, 0.0060001800011014146423}, {-0.005,0.0050001041671093232710},
{-0.004, 0.0040000533334783891518}, {-0.003, 0.0030000225000344235520},
{-0.002, 0.0020000066666711999153}, {-0.001, 0.0010000008333334749993},
{0., 0},
{0.001, -0.0010000008333334749993}, {0.002, -0.0020000066666711999153},
{0.003, -0.0030000225000344235520}, {0.004, -0.0040000533334783891518}}  *)

ListPlot[soln]


Here are the solutions for 50 zeros

fun = Sinh[y] - (y + Tanh[x]) Cosh[y] - y Sech[x]^2;

xlist = Array[# &, 50, {-5, 5}];

sol = NSolve[fun /. x :> #, {y}, Reals][[1, 1, 2]] & /@ xlist;

zero = Transpose[{sol, Array[0 &, Length@sol]}];

Plot[Evaluate@Table[fun, {x, {xList}}], {y, -2.5, 2.5},
Epilog -> {Red, PointSize@Large, Point@zero},
ImageSize -> Large,
PlotRange -> {Automatic, {-2, 2}},
PlotStyle -> Gray]