Order of error in reusing numerical derivative?

When using numerical approximations of the taylor series form

$$f_N(x \pm h)=\sum_{k=0}^{N}{\frac{(\pm h)^kf^{(k)}(x)}{k!}}$$

for derivatives with error of order $O(h^2)$

$$y'=\frac{y(x+h)-y(x-h)}{2h}$$ and

$$y''=\frac{y(x+h)-2y(x)+y(x-h)}{h^2}$$

How would the error change (if it does) if you were to reuse this equation to approximate another derivative?

Specifically, i have this problem $y'''=(y'')^2$, which can be simplified into $u'=u^2$ where $u=y''$. For notation $u(x\pm h)=u_{\pm h}$, i've simplified this equation into $\frac{u_{h}-u_{-h}}{2h}=u_{0}^2$ which i assume holds the error of $O(h^2)$, but then i'm not sure if further expanding this into $u_{0}=\frac{y_h-2y_0+y_{-h}}{h^2}$ will net the same error? I've simulated this for a specific intial condition case in matlab and the solution for $y$ seems to become unstable very quickly.

• What are the initial conditions you use for $u$ (and $y$) and for what $t$ do you see an instabillity? What is the stepsize you use? The analytial solution to $u' = u^2$ is $u(t) = \frac{u(0)}{1 - u(0) t}$ which has a singularity at $t = \frac{1}{u_0}$ so this could be some of the reason. – Winther Oct 5 '16 at 17:22
• I'm using $y(0)=y'(0)=0, y''(0)=-1$ which gives me $y(x)=x-(x+1)ln(x+1), y'(x)=-ln(x+1), y''(x)=\frac{-1}{x+1}$. I'm using the exact values from $y$ to put into the iterative solution for $y_h$ for $x \geq 0$. I'm not sure if it is possible to apply the approximation twice, but the iterative solution i ended up getting was $y_h=2(y_0-y_{-2h})+y_{-3h}+\frac{2}{h}(y^2_0+y^2_{-2h}+4(y^2_{-h}-y_{-h}(y_0+y_{-2h})-y_0y_{-2h}))$, but i'm not sure if it is allowable to use the approximation twice and then shift the index, unless the problem is the initial condition – Baklava Gain Oct 7 '16 at 13:06
• Trying to discretize $y''' = (y'')^2$ leads to a very complicated expression. I would instead solve $u' = u^2$ and then solve $y'' = u$. This can in fact be solved for simulatanously by defining the 3D vector $\vec{Y}$ with $Y_1 = y''$, $Y_2 = y'$ and $Y_3 = y$ then the system can be written on the form $\frac{d}{dx}\vec{Y} = \vec{f}(\vec{Y})$ as $\frac{d}{dx}\pmatrix{Y_1\\Y_2\\Y_3} = \pmatrix{Y_1^2\\Y_1\\Y_2}$. I tried solving the system with Euler's method with $n=100$ points between $x=0$ and $x=1$ and it worked perfectly. – Winther Oct 7 '16 at 15:50
• Ok, i got the solution after using this with the euler method and this method as well.. I'm assuming that the system of equations $y_1'=y^2, y_2'=y_1, y_3'=y_2$ cannot be combined into a single equation which i think was my initial problem, each equation has to be solved separately. Thanks for the help. – Baklava Gain Oct 8 '16 at 2:29