Having problem with empty matrix...Can anyone help?

here is the code..

   disp('dy/dx=1+xy');
h=0.2;
xfinal=0.6;
y(1)=2;
x(1)=0;
f=@(x,y) 1+x*y;
for i=1:ceil(xfinal/h)
    x(i+1)=x(i)+h;
    k1=f((x(i)),y(i));
    k2=f((x(i)+.5*h),(y(i)+.5*h*k1));
    k3=f((x(i)+.5*h),(y(i)+.5*h*k2));
    k4=f((x(i)+h),(y(i)+k3*h));
    y(i+1)=y(i)+((1/6)*(k1+(2*k2)+(2*k3)+k4)*h);
end
plot(x,y);
xlabel('x');
ylabel('y');
grid on;
idx1=find(x==0.2);
y1=y(idx1)
idx2=find(x==0.4);
y2=y(idx2)
idx3=find((x==0.6))
y3=y(idx3)disp('dy/dx=1+xy');

it returns an empty matrix error for y3...what should i do?

Answer 1

[Star Strider]

The problem is that you are encountering floating point approximation error.

If you do the comparison:

x_error = x(4) - 0.6

x_error =
   111.0223e-018

The solution is to allow for the tolerance. The easiest way is to change your ‘idx3’ assignment to:

idx3=find((x<=0.6),1,'last')

which works.

For a thorough discussion, see Why is 0.3 - 0.2 - 0.1 (or similar) not equal to zero?

Answer 2

[Walter Roberson]

Do not compare floating point numbers for exact equality. Compare them with a tolerance instead. Two different ways of calculating the same result can give different results with floating point. For example 1 + 1E200 - 1E200 will not give 1