Problem with runge-kutta adaptive algorithm

Hey,

I have the following problem. I'm trying to write a program in Matlab, that would implement Runge-Kutta 2 algorithm, but with changing step size, so the adaptive one. In the main script I need to write a code based on this pseudocode:

    while(t<tk)
    h=h_start; // h-step size

    while(crit>eps)
    y_next=RK2(t, y, h); //here RK2 method is called
    y_next_2=RK2(t, y, h/2);
    y_next_2=RK2(t, y_next_2, h/2);

    crit=norm(y_next-y_next_2);

    if(crit>eps)
    h=h/2;
    end
  y(t+1)=y_next_2;
 end
 end

So, basicaly what's need to be done, is to draw one point of a diagram in each loop, changing the step size h until criterium (crit) is less than epsilon (eps).

Below, there is the code I've already done in Matlab. The problem is, in this code I'm not getting one diagram but few of them, because the plot function is in the loop of the RK2 method. I couldn't find any other way to send the values of the function to the main script.

Here's the script:

Main:

   h_start=5;
   %h_start=0.1;
   h=h_start;

   crit=1;
   eps=0.00001;
   q=1;

while(crit>eps)
 y_next=RK2(h);
 y_next_2=RK2(h/2);

   crit=norm(y_next_2-y_next);

   if(crit>eps)
   h=h/2;
   end

   q=q+1;

 end

%plot(t, y2, 'r');
q

RK2:

   function [y1, t] = RK2(h)

   y1(1)=0; %initial conditions
   y2(1)=0;
   p(1)=0;

   c1=0.26;
   c2=0.1;
   c3=1;

   a=0.13;
   b=0.013;

   F_xy= @(t, y1, y2, p) c1*y1*(y1-a)*(1-y1)-c2*y2+p;
   G_xy= @(t, y1, y2) b*(y1-c3*y2);

   t(1)=0;

   n=80;
   %n=2000;

   for j = 1:n
   if(t(j)>50 && t(j)<60) 
    p(j)=0.05;
   else
    p(j)=0;
   end

   k1 = h * F_xy(t(j), y1(j), y2(j), p(j));
   l1= h * G_xy(t(j), y1(j), y2(j));
   k2 = h * F_xy( t(j) + h/2, y1(j) + k1/2, y2(j) + k1/2, p(j));
   l2 = h * G_xy( t(j) + h/2, y1(j) + l1/2, y2(j) + l1/2);

   y1(j+1) = y1(j) + (k1 + k2)/2;
   y2(j+1) = y2(j) + (l1 + l2)/2;

   t(j+1)=t(j)+h;
   hold on
   plot(t, y1, 'r');

 end

 end

[Jan]

Please explain, why you enclose the code blocks in "if true". It is a good idea to omit everything, which is not connected to the problem. Adding this useless line even to the pseudo-code version looks strange.

Usually "adaptive RK integrators" adjust the step size in each step. Your code is not "adaptive" in the usual sense.

Your pseudo-code is hard to understand. "while(t<tk)" requires that "t" is changes inside the loop, but it isn't in your code. "y_nast_2" is not defined. "while(crit>eps)" will suffer from an undefined "crit" in the first iteration. Even it is only pseudo-code, it must be correct also. Otherwise it will not help to create a correct program.

[Daniel]

Yes, you're right it's useless. I used "code" button when I wrote the message- I thought it would put the code in a distinctive part.

[Daniel]

I corrected "y_nast_2", it should be "y_next_2". And "t" is changing in this line: "t(j+1)=t(j)+h". This is the main problem that I had. I don't know how to include t in the main program, so that the rk2 method would still work. "crit" variable I've defined in main program. In pseudocode it's not defined.

I used t variable in rk2, because it's used to estimate the next value of coefficients k1 and k2 and then to estimate value of "y1(j+1)". The exact lines, where its used:

k1 = h * F_xy(t(j), y1(j), y2(j), p(j));

k2 = h * F_xy( t(j) + h/2, y1(j) + k1/2, y2(j) + k1/2, p(j));