# A code for producing system of equations for central difference method

I am trying to solve the differential equation $$y''(x) = y(x)+\sin[y'(x)]$$ numerically over the interval $$[0,1].$$ Now using central difference method I arrived at the equations of the type $$y_k = \frac{1}{h^2+2} \left(y_{k-1} + y_{k-1} - h^2\sin\left(\frac{y_{k+1}-y_{k-1}}{2h}\right)\right)$$

for $$k=1,2,\ldots, m-1.$$ I am trying to solve for $$m=10$$ making a system of equations with $$y_0=0$$ and $$y_{10}=1.$$ The issue is I am not able to produce the equations in Mathematica.

SysEq[xs_List] := (1/(h^2 + 2) (#[[1]] + #[[3]] - h^2*Sin[(#[[3]] - #[[1]])/(2*h)])) & /@ Partition[xs, 10, 1, {1, 1}];


So that I can call it

SysEq[{x1, x2, x3, x4, x5, x6, x7, x8, x9, x10}]


or may be an easy way to create whole system of equations.

to generate the equations, just use Table command. To solve them is different issue. Since they are non-linear.

$$y^{\prime\prime}=y+\sin\left( y^{\prime}\right)$$

3 point centered difference for $$y^{\prime\prime}$$ is

$$y^{\prime\prime}=\frac{1}{h^{2}}\left( y_{n-1}-2y_{n}+y_{n+1}\right)$$

centered difference for $$y^{\prime}$$ is

$$y^{\prime}=\frac{1}{2h}\left( y_{n+1}-y_{n-1}\right)$$

The ode becomes

$$\frac{1}{h^{2}}\left( y_{n-1}-2y_{n}+y_{n+1}\right) =y_{n}+\sin\left( \frac{1}{2h}\left( y_{n+1}-y_{n-1}\right) \right)$$

let $$L=1$$ and 11 points. 9 are internal points, and 2 are boundary points. Hence number of intervals is 10.

Therefore $$h=\frac{1}{10}$$ and $$n=0,1,\cdots,10$$. Where $$y_{0}=0,y_{10}=1$$. For $$n=1$$

$$\frac{1}{h^{2}}\left( y_{0}-2y_{1}+y_{2}\right) =y_{1}+\sin\left( \frac {1}{2h}\left( y_{2}-y_{0}\right) \right)$$

For $$n=2$$

$$\frac{1}{h^{2}}\left( y_{1}-2y_{2}+y_{3}\right) =y_{2}+\sin\left( \frac {1}{2h}\left( y_{3}-y_{1}\right) \right)$$

And so on. For $$n=9$$

$$\frac{1}{h^{2}}\left( y_{8}-2y_{9}+y_{10}\right) =y_{9}+\sin\left( \frac {1}{2h}\left( y_{10}-y_{8}\right) \right)$$

Hence the equations are

\begin{align*} y_{0} & =0\\ \frac{1}{h^{2}}\left( y_{0}-2y_{1}+y_{2}\right) & =y_{1}+\sin\left( \frac{1}{2h}\left( y_{2}-y_{0}\right) \right) \\ \frac{1}{h^{2}}\left( y_{1}-2y_{2}+y_{3}\right) & =y_{2}+\sin\left( \frac{1}{2h}\left( y_{3}-y_{1}\right) \right) \\ & \vdots\\ \frac{1}{h^{2}}\left( y_{8}-2y_{9}+y_{10}\right) & =y_{9}+\sin\left( \frac{1}{2h}\left( y_{10}-y_{8}\right) \right) \\ y_{10} & =1 \end{align*}

Code is

numberInternalPoints = 9;
h                    = 1/(numberInternalPoints+1);
y[0]                 = 0;
y[numberInternalPoints + 1] = 1;

vars = Table[y[n], {n, 1, numberInternalPoints}]
eqs = Table[
1/h^2*(y[n - 1] - 2 y[n] + y[n + 1]) ==
y[n] + Sin[1/(2*h)*(y[n + 1] - y[n - 1])], {n,1, numberInternalPoints}];
TableForm[eqs]


9 equations, 9 unknowns.

Update

As suggested by user64494 below, FindRoot can be used to solve these equations, and the result agree with NDSolve

sol = FindRoot[
eqs, {{y[1], 0}, {y[2], 0}, {y[3], 0}, {y[4], 0}, {y[5], 0}, {y[6],
0}, {y[7], 0}, {y[8], 0}, {y[9], 0}}]
data = Table[{n*h, Last@sol[[n]]}, {n, 1, numberInternalPoints}]
data = PrependTo[data, {0, 0}];
data = AppendTo[data, {1, 1}]


Gives

{{0, 0}, {1/10, 0.0558825}, {1/5, 0.117883}, {3/10, 0.187164}, {2/5,
0.265029}, {1/2, 0.352914}, {3/5, 0.452386}, {7/10, 0.565109}, {4/5,
0.692812}, {9/10, 0.837223}, {1, 1}}


And

p = ListLinePlot[data, Mesh -> All, MeshStyle -> Red]
ClearAll[y, x]
sol = NDSolveValue[{y''[x] == y[x] + Sin[y'[x]], y[0] == 0,
y[1] == 1}, y, {x, 0, 1}];
p2 = Plot[sol[x], {x, 0, 1}, PlotStyle -> Blue];
Show[p, p2]


• A spot on the Sun: NSolve[eqs, vars, Reals] generates "NSolve::nsmet: This system cannot be solved with the methods available to NSolve". Commented Sep 21, 2023 at 19:21
• Just to compare. The command of Maple DirectSearch:-SolveEquations([100*(-2*y[1]+y[2]) = sin(5*y[2])+y[1], 100*(y[1]-2*y[2]+y[3]) = sin(5*(-y[1]+y[3]))+y[2], 100*(y[2]-2*y[3]+y[4]) = sin(5*(-y[2]+y[4]))+y[3], 100*(y[3]-2*y[4]+y[5]) = sin(5*(-y[3]+y[5]))+y[4], 100*(y[4]-2*y[5]+y[6]) = sin(5*(-y[4]+y[6]))+y[5], 100*(y[5]-2*y[6]+y[7]) = sin(5*(-y[5]+y[7]))+y[6], 100*(y[6]-2*y[7]+y[8]) = sin(5*(-y[6]+y[8]))+y[7], 100*(y[7]-2*y[8]+y[9]) = sin(5*(-y[7]+y[9]))+y[8], 100*(1+y[8]-2*y[9]) = sin(5*(1-y[8]))+y[9]]) results in Commented Sep 21, 2023 at 19:46
• [y[1] = 0.558824815612752e-1, y[2] = .117882524237100, y[3] = .187164143110649, y[4] = .265028696216958, y[5] = .352914417317458, y[6] = .452385864615537, y[7] = .565109488648346, y[8] = .692812302307770, y[9] = .837223071943341] in a moment. Commented Sep 21, 2023 at 19:47
• FindRoot[{100 (-2 y[1] + y[2]) == Sin[5 y[2]] + y[1], 100 (y[1] - 2 y[2] + y[3]) == Sin[5 (-y[1] + y[3])] + y[2], 100 (y[2] - 2 y[3] + y[4]) == Sin[5 (-y[2] + y[4])] + y[3], 100 (y[3] - 2 y[4] + y[5]) == Sin[5 (-y[3] + y[5])] + y[4], 100 (y[4] - 2 y[5] + y[6]) == Sin[5 (-y[4] + y[6])] + y[5], 100 (y[5] - 2 y[6] + y[7]) == Sin[5 (-y[5] + y[7])] + y[6], 100 (y[6] - 2 y[7] + y[8]) == Sin[5 (-y[6] + y[8])] + y[7], 100 (y[7] - 2 y[8] + y[9]) == Sin[5 (-y[7] + y[9])] + y[8], 100 (1 + y[8] - 2 y[9]) == Sin[5 (1 - y[8])] + y[9]},  Commented Sep 21, 2023 at 20:03
• {{y[1], 0}, {y[2], 0}, {y[3], 0}, {y[4], 0}, {y[5], 0}, {y[6], 0}, {y[7], 0}, {y[8], 0}, {y[9], 0}}] outputs {y[1] -> 0.0558825, y[2] -> 0.117883, y[3] -> 0.187164, y[4] -> 0.265029, y[5] -> 0.352914, y[6] -> 0.452386, y[7] -> 0.565109, y[8] -> 0.692812, y[9] -> 0.837223} Commented Sep 21, 2023 at 20:04