Wrong value for the derivative of the solution of an ODE in MATLAB

Question

Trying to solve a Cauchy problem for a nonlinear ODE with the classical Runge-Kutta method described in here (see, ode4). After running this part

q = 10;
T = 5;
N = @(g1, g2) g1^2 + sin(g2);

t = linspace(0, T, q);
GCC = nonlinearGreenCC(N, 1, T);
di = diff(GCC);

I evaluate GCC(1) and di(1). While GCC(1) = 0 as expected, di(1) = 1.6e-05. Can't understand why, because the Cauchy condition on the first derivative is 1. How to fix the inaccuracy/mistake? Any help is highly appreciated.

The function nonlinearGreenCC is as follows:

function G = nonlinearGreenCC(N, a0, T)

h = .0001;
eqGreen = @(t, g)[g(2); - N(g(1), g(2))];
Cc = [0, a0];
G = ode4(eqGreen, 0, h, T, Cc);

end

Answer 1

GCC is an array of pairs of values, starting for a more exact integration with something like

[  0.00000000e+00,   1.00000000e+00],
[  9.99957927e-05,   9.99915855e-01],
[  1.99983171e-04,   9.99831715e-01],
[  2.99962136e-04,   9.99747579e-01],
[  3.99932687e-04,   9.99663448e-01],

for exact values, but your result using the ode4 function starts with something like

 0                         1
 1.66652642851402e-05      1.00001666526429
-1.40231141152723e-05      0.999985976885885
 1.66638620438405e-05      1.00001666386204
-1.40217119017434e-05      0.999985978288098

It does not contain the sampling of the time interval, indeed on the level that you call diff, the value of h is unknown. There is no way that diff can compute the difference quotients that you seem to expect. If you read the documentation you should find that it computes just the differences. And the first difference returns just the first value 1.66652642851402e-05.

---

The direct error that causes the ode4 algorithm to produce strange results is that eqGreen produces column vectors where it should return row vectors. As the initial value is a row vector, the result of adding row and column vectors in ode4 produces a 2x2 matrix that is added as 2 rows to the result with confusing consequences. Having both row vectors puts the results into one row with alternating value and derivative. Change to

eqGreen = @(t, g)[g(2), - N(g(1), g(2))];
Cc = [0, a0];

and the result is

GCC(1:5,:) =

                     0                     1
  9.99957927208439e-05     0.999915855174488
  0.000199983171186392     0.999831714893813
  0.000299962135851046     0.999747579156327
  0.000399932687169042     0.999663447960381

di(1:5,:)=

  9.99957927208439e-05  -8.41448255122224e-05
  9.99873784655483e-05  -8.41402806746050e-05
  9.99789646646540e-05  -8.41357374861129e-05
  9.99705513179961e-05  -8.41311959463020e-05
  9.99621384254097e-05  -8.41266560547282e-05

If you scale the last by dividing by h=1e-4 you get the expected results.