Trying to change the output vector for Heun's non-self-starting

I have created a program to run for non-self-starting Heun's method for a given problem in class. But what I am noticing is that my "y" output has one extra 0 at the start that id like to remove.

This is the given problem

And this is my code

f= @(t,y)(sin(y)/y+2*t);       
a = input('Enter LOWER limit of domain, a=  ');
b = input('Enter UPPER limit of domain, b=  ');
h = input('Enter step size, h= ');               
t=a:h:b;
y=zeros(size(t));
y(2) = input('Enter Initial condition of y(0),=  ');
y(1) = input('Enter Initial condition of y(i-1),=  ');
n=numel(t);
for i=2:n
nil(i-1)=y(i-1)+(2*h*f(t(i-1),y(i)));
y(i+1)=y(i)+(h/2)*(f(t(i),nil(i-1))+f(t(i-1),y(i)));
end

Output is as follows (im only interested in the "t" andn "y") Now I very well could be wrong in this but given the initial condition of the problem, should my first "y" output be 1 and not 0? and if so, how can i fix that.

[Johan]

The first two values of y are your initial conditions:

y(2) = input('Enter Initial condition of y(0),=  ');
y(1) = input('Enter Initial condition of y(i-1),=  ');

According to your problem they should be 0 and 1 if I understand correctly.

I'm not familiar with non self starting Heun's method but the solution may be to shift your vector t by -1 step ?

t=(a-h):h:b;

Answer 1

[Shivam]

Hi Paul,

As per my understanding, you are having an issue filling output vector 'y' while solving given differential equation with non-self-starting Heun's method.

You can follow the below workaround to resolve the issue:

• Modify vector t to have its first value as -1.

     t = [-1, a:h:b];

• Modify for loop to calculate the correct predictor and corrector following the non-self-starting Heun's method.

    % ..
    % ...
    n=numel(t);
    for i=2:n-1
        nil = y(i-1)+f(t(i), y(i))*2*h;                     % Predictor
        y(i+1) = y(i)+(h/2)*(f(t(i), y(i))+f(t(i+1), nil)); % Corrector
    end

The first value of vector y corresponds to y(t(1)) ~= y(-1) = 0.

In case you need the output vector y with no output value corresponding to t = -1, you can modify y after above calculations:

    y = y(:, 2:end)

I hope it helps in resolving the issue.