RK 4 method for ODES

When I run this code it gives me error. I tried to change some value but no luck

function RK4(a,b,N)

if nargin<3, a=0; b=1; N=20; end

h=(b-a)/N;

y(1)=1;

y1(1)=2;

for i=1:N

t=(i-1)*h;

k1=f(t,y(i),y1(i));

k2=f(t+h/2,y(i)+k1*h/2,y1(i)+k1*h/2);

k3=f(t+h/2,y(i)+k2*h/2,y1(i)+k2*h/2);

k4=f(t+h,y(i)+k3*h,y1(i)+k3*h);

y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4)*h;

end

for i=1:N+1

t=(i-1)*h;

exact(i)=exp(t)+t;

end

error=max(abs(exact-y))

disp(y(N));

figure;

plot(i,y,'LineWidth',2);

xlabel('X');

ylabel('Y');

end

function fty=f(t,y,y1)

fty=-y1+exp(-7*t)*(t-y)+2*exp(t)+exp(-6*t)+1;

end

Output

Index exceeds the number of array elements. Index must not exceed 1.

Error in RK4 (line 28)

k1=f(t,y(i),y1(i));

Answer 1

[Steven Lord]

y(1)=1;
y1(1)=2;
for i=1:N

At the start of the loop (during iteration 1) y and y1 each have 1 element.

    t=(i-1)*h;
    k1=f(t,y(i),y1(i));

At iteration 1 this line uses element 1 of the y and y1 variables.

At iteration 2 this line uses element 2 of the y and y1 variables.

    k2=f(t+h/2,y(i)+k1*h/2,y1(i)+k1*h/2);
    k3=f(t+h/2,y(i)+k2*h/2,y1(i)+k2*h/2);
    k4=f(t+h,y(i)+k3*h,y1(i)+k3*h);
    y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4)*h;

At iteration 1 this line assigns to element 2 of y, so at the start of iteration 2 y will have 2 elements.

Nowhere in this loop do you assign to any element of y1, so it will always remain 1 element long. This is the cause of the error at iteration 2. It doesn't make sense to ask for the value of the second element of a vector with one element.

end

If your goal is to solve a system of ODEs, use one of the ones included in MATLAB rather than writing your own.

If your goal is to implement your own ODE solver (because it's a homework assignment, for example) you need to determine what values to assign to subsequent elements in y1 inside the loop like you're assigning to y.

[Syed Abdul Rafay]

y1’ = - y1(t) + e-7t (t - y) + 2 et + e-6t + 1

This is my equation

[Syed Abdul Rafay]

I have three unknow and I have two conditions y(0) = 0 and y1(0) = 2 the y(t) and y1(t). So I want to solve for y1 also

[Syed Abdul Rafay]

in the first calculation the loop should take value of that of y(1) and y1(1) the y(2) and y(2) value will be gained from the end of the loop of first iteration.

[Steven Lord]

the y(2) and y(2) value will be gained from the end of the loop of first iteration.

You assign to y(2) in your loop [technically y(i+1) where i is equal to 1.] You never assign to y1(2).

[Syed Abdul Rafay]

I corrected my code thanks for the help.