# A way to generalize my code for this non linear ODE problem

I am trying to solve the differential equation $$y''(x) = y(x)+\sin[y'(x)]$$ using "Fixed Point Iteration" over the interval $$[0,1].$$ Now using central difference method I arrived at the equations of the type $$F(y_0,y_1,\ldots,y_k) : 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) =0$$

for $$k=1,2,\ldots, m-1.$$ I am trying to solve for $$m=50$$ making a system of equations with $$y_0=0$$ and $$y_{50}=1.$$ Then convert it into fixed point equations $$G(y_0,y_1,\ldots,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) = y_k$$

The issue is I am not able to produce the equations (A general way) in Mathematica. I have this code for another problem. The first two lines are the one I need udated for this example, so that they produce the equired $$F's$$ and $$G's$$ for any $$m.$$

    Quit
FList[xs_List] := (#[[1]]^2 #[[2]] - 1) & /@ Partition[xs, 2, 1, {1, 1}];
GList[xs_List] := RotateLeft[1/Sqrt[xs]];
JacMat[xs : {Symbol__}] := D[#, {xs}] & /@ GList[xs];

Module[{eps = 10^-160, max = 15, dim = 5, per = 240,  a, b, k, ts, xs,
x}, ts = SetAccuracy[ConstantArray[11/10, dim], per];
xs = Array[x, dim];
a = 1;
b = Norm[FList[ts], Infinity];
k = 1;
While[k <= max && a + b > eps,
M = SetAccuracy[-JacMat[xs] /. Thread[xs -> ts], per];
J = SetAccuracy[
Inverse[IdentityMatrix[dim] - JacMat[xs] /. Thread[xs -> ts]],
per];
hs = GList[GList[J . M . ts + J . GList[ts]]];
a = SetAccuracy[Norm[(ts - hs), Infinity], per];
b = SetAccuracy[Norm[FList[hs], Infinity], per];
Print["k=", NumberForm[k, 20]];
Print["hs=", NumberForm[hs, 20]];
Print["a=", NumberForm[a, 20]];
Print["b=", NumberForm[b, 20]];
ts = hs;
k = k + 1
];
Print["The number of iterations is:" NumberForm[(k - 1), 2]];
Print["a+b=" NumberForm[a + b, 20]]]

• Is this a typo that you have: y[k-1]+y[k-1] in F and G? Commented Sep 23, 2023 at 16:52
• yes thats a typo I have fixed it. Thanks for pointing it out Commented Sep 23, 2023 at 17:01
• I do not understand what the difference between F and G is. It is the same equation. Commented Sep 23, 2023 at 17:36
• If I explain in one variable The F(x) = 0 is conveted into G(x) = x,(Fixed point problem) in this case just by rearrangement. Commented Sep 23, 2023 at 17:42

Here is an example how to solve for F.

f[ik_] :=  ylist[k] -  1/(h^2 + 2) (ylist[k + 1] + ylist[k - 1] -
h^2 Sin[(ylist[k + 1] - ylist[k - 1])/(2 h)]) == 0

Then for a numerical example, using FindInstance and e.g. n=10 and h=1:

n = 10; h = 1;
eq = Table[f[k], {k, 1, n - 1}] /. {ylist[0] -> 0, ylist[n] -> 1};
vs = Table[ylist[i], {i, 1, n}];
FindInstance[eq, vs[[1 ;; n - 1]]] // N

{{ylist[1.] -> 2.20295*10^-7, ylist[2.] -> 1.32177*10^-6,
ylist[3.] -> 7.26974*10^-6, ylist[4.] -> 0.0000396531,
ylist[5.] -> 0.000216109, ylist[6.] -> 0.0011777,
ylist[7.] -> 0.00641784, ylist[8.] -> 0.0349724,
ylist[9.] -> 0.190322}}

• Thank you for the answer. But thats not what I am looking for. I have my own ways to solve the equations. I just need a way to produce the equations (first three lines). so that I can run them through my codes (second part of the code above) Commented Sep 24, 2023 at 4:39
• @Learner Look for NDSolveFiniteDifferenceDerivative` to get the discretization. Commented Sep 24, 2023 at 10:26
• @UlrichNeumann I can not use that. Commented Sep 24, 2023 at 12:25
• @Learner Why not? Commented Sep 24, 2023 at 14:31
• I have my own ways to solve the equations. I just need a way to produce the equations (first three lines). so that I can run them through my codes (second part of the code above) Commented Sep 24, 2023 at 17:03