Matlab Plot can't work.

Hi Guys,

Below is my code, and I have a quesion, when I modify the plot function from plot(x,z,'-*r') to plot(x,z), the figure can't show, Does anyone know why?

clc

clear

clc

f = @(x,y,z) 1/y*(y+z); %where k = 1 for y''

g = @(x,y,z) z ; %for y'

h = 0.1;

x= 0;

z= 1;

y= 1;

%value of y at x = %0.2f and z = 0.2%f is %f\n

for i = 1:25

k1 = h.*f(x,y,z);

j1 = h.*g(x,y,z);

k2 = h.*f(x+h/2,y+k1/2,z+j1/2);

j2 = h.*g(x+h/2,y+k1/2,z+j1/2);

k3 = h.*f(x+h/2,y+k2/2,z+j2/2);

j3 = h.*g(x+h/2,y+k2/2,z+j2/2);

k4 = h.*f(x+h,y+k3,z+j3);

j4 = h.*g(x+h,y+k3,z+j3);

k = (k1+2*k2+2*k3+k4)./6;

j = (j1+2*j2+2*j3+j4)./6;

x = x + h;

z = z + j;

y = y + k;

fprintf('i = %8.0f, y = %8.6f, x = %8.2f, z = %8.6f\n',i,y,x,z);

plot(x,z,'-*r')

hold on

plot(x,y,'-ob')

end

Answer 1

[KSSV]

clc
clear
clc
f = @(x,y,z) 1/y*(y+z);  %where k = 1   for y''
g = @(x,y,z) z ;   %for y'


h    = 0.1;
x= 0;
z= 1;
y= 1;
%value of y at x = %0.2f and z = 0.2%f is %f\n
figure
hold on
for i = 1:25
    
    k1 = h.*f(x,y,z);
    j1 = h.*g(x,y,z);
    k2 = h.*f(x+h/2,y+k1/2,z+j1/2);
    j2 = h.*g(x+h/2,y+k1/2,z+j1/2);
    k3 = h.*f(x+h/2,y+k2/2,z+j2/2);
    j3 = h.*g(x+h/2,y+k2/2,z+j2/2);
    k4 = h.*f(x+h,y+k3,z+j3);
    j4 = h.*g(x+h,y+k3,z+j3);
    k  = (k1+2*k2+2*k3+k4)./6;
    j  = (j1+2*j2+2*j3+j4)./6;
    x = x + h;
    z = z + j;
    y = y + k;
    fprintf('i = %8.0f, y = %8.6f, x = %8.2f, z = %8.6f\n',i,y,x,z);
    
    plot(x,z,'-*r')    
    plot(x,y,'-ob')
end

i =        1, y = 1.195848, x =     0.10, z = 1.105171
i =        2, y = 1.385917, x =     0.20, z = 1.221403
i =        3, y = 1.572763, x =     0.30, z = 1.349858
i =        4, y = 1.758013, x =     0.40, z = 1.491824
i =        5, y = 1.942803, x =     0.50, z = 1.648721
i =        6, y = 2.127991, x =     0.60, z = 1.822118
i =        7, y = 2.314263, x =     0.70, z = 2.013752
i =        8, y = 2.502200, x =     0.80, z = 2.225540
i =        9, y = 2.692311, x =     0.90, z = 2.459601
i =       10, y = 2.885062, x =     1.00, z = 2.718280
i =       11, y = 3.080891, x =     1.10, z = 3.004163
i =       12, y = 3.280220, x =     1.20, z = 3.320114
i =       13, y = 3.483461, x =     1.30, z = 3.669293
i =       14, y = 3.691027, x =     1.40, z = 4.055196
i =       15, y = 3.903333, x =     1.50, z = 4.481684
i =       16, y = 4.120803, x =     1.60, z = 4.953026
i =       17, y = 4.343871, x =     1.70, z = 5.473940
i =       18, y = 4.572985, x =     1.80, z = 6.049639
i =       19, y = 4.808610, x =     1.90, z = 6.685885
i =       20, y = 5.051229, x =     2.00, z = 7.389045
i =       21, y = 5.301345, x =     2.10, z = 8.166157
i =       22, y = 5.559487, x =     2.20, z = 9.024998
i =       23, y = 5.826203, x =     2.30, z = 9.974165
i =       24, y = 6.102073, x =     2.40, z = 11.023156
i =       25, y = 6.387702, x =     2.50, z = 12.182471

The suggested way is this:

clc
clear
clc
f = @(x,y,z) 1/y*(y+z);  %where k = 1   for y''
g = @(x,y,z) z ;   %for y'


h    = 0.1;
N = 25 ;
x = zeros(N,1) ;
y = zeros(N,1) ;
z = zeros(N,1) ;

x(1)= 0;
z(1)= 1;
y(1)= 1;
%value of y at x = %0.2f and z = 0.2%f is %f\n
for i = 2:25
    
    k1 = h*f(x(i-1),y(i-1),z(i-1));
    j1 = h*g(x(i-1),y(i-1),z(i-1));
    k2 = h*f(x(i-1)+h/2,y(i-1)+k1/2,z(i-1)+j1/2);
    j2 = h*g(x(i-1)+h/2,y(i-1)+k1/2,z(i-1)+j1/2);
    k3 = h*f(x(i-1)+h/2,y(i-1)+k2/2,z(i-1)+j2/2);
    j3 = h*g(x(i-1)+h/2,y(i-1)+k2/2,z(i-1)+j2/2);
    k4 = h*f(x(i-1)+h,y(i-1)+k3,z(i-1)+j3);
    j4 = h.*g(x(i-1)+h,y(i-1)+k3,z(i-1)+j3);
    k  = (k1+2*k2+2*k3+k4)./6;
    j  = (j1+2*j2+2*j3+j4)./6;
    x(i) = x(i-1) + h;
    z(i) = z(i-1) + j;
    y(i) = y(i-1) + k;
    fprintf('i = %8.0f, y = %8.6f, x = %8.2f, z = %8.6f\n',i,y(i),x(i),z(i));
end

i =        2, y = 1.195848, x =     0.10, z = 1.105171
i =        3, y = 1.385917, x =     0.20, z = 1.221403
i =        4, y = 1.572763, x =     0.30, z = 1.349858
i =        5, y = 1.758013, x =     0.40, z = 1.491824
i =        6, y = 1.942803, x =     0.50, z = 1.648721
i =        7, y = 2.127991, x =     0.60, z = 1.822118
i =        8, y = 2.314263, x =     0.70, z = 2.013752
i =        9, y = 2.502200, x =     0.80, z = 2.225540
i =       10, y = 2.692311, x =     0.90, z = 2.459601
i =       11, y = 2.885062, x =     1.00, z = 2.718280
i =       12, y = 3.080891, x =     1.10, z = 3.004163
i =       13, y = 3.280220, x =     1.20, z = 3.320114
i =       14, y = 3.483461, x =     1.30, z = 3.669293
i =       15, y = 3.691027, x =     1.40, z = 4.055196
i =       16, y = 3.903333, x =     1.50, z = 4.481684
i =       17, y = 4.120803, x =     1.60, z = 4.953026
i =       18, y = 4.343871, x =     1.70, z = 5.473940
i =       19, y = 4.572985, x =     1.80, z = 6.049639
i =       20, y = 4.808610, x =     1.90, z = 6.685885
i =       21, y = 5.051229, x =     2.00, z = 7.389045
i =       22, y = 5.301345, x =     2.10, z = 8.166157
i =       23, y = 5.559487, x =     2.20, z = 9.024998
i =       24, y = 5.826203, x =     2.30, z = 9.974165
i =       25, y = 6.102073, x =     2.40, z = 11.023156

plot(x,z,'-*r')
hold on
plot(x,y,'-ob')

[Jason Ciou]

Hi @KSSV,

May I know why I need to start at 2?

[KSSV]

You are using previous value right?

[Jason Ciou]

Yes