clc;
clear;
xdel(winsid());

// Initializing variables
h = 0.05;
x = 3;
y = 1;
z = 2;
xt = 0.2;
yt = 0.1;
zt = 0.05;

// Start values equals zero
for a = 1:2000
 x_st = 0;
 y_st = 0;
 z_st = 0;
 xt_st = 0;
 yt_st = 0;
 zt_st = 0;
end

// myval equals zero in all dimensions
for i = 1:2000
 for j = 1:3
    myvalue = 0;
 end;
end

// Using RK4 method to solve equations for 2000 points
for i = 1:2000
 myvalue(i,1) = y + z;
 myvalue(i,2) = -x + 0.5*y;
 myvalue(i,3) = x^2 - z;

 kx1 = h*(y + z);
 ky1 = h*(-x + 0.5*y);
 kz1 = h*(x^2 - z);

 kx2 = h*((y+ky1)/2 + (z+kz1)/2);
 ky2 = h*((-x+kx1)/2 + 0.5*((y+ky1)/2));
 kz2 = h*(((x+kx1)^2)/2 - (z+kz1)/2);

 kx3 = h*((y+ky2)/2 + (z+kz2)/2);
 ky3 = h*((-x+kx2)/2 + 0.5*((y+ky2)/2));
 kz3 = h*(((x+kx2)^2)/2 - (z+kz2)/2);

 kx4 = h*((y+ky3) + (z+kz3));
 ky4 = h*((-x+kx3) + 0.5*(y+ky3));
 kz4 = h*((x+kx3)^2 - (z+kz3));

 xt = (kx1+2*kx2+2*kx3+kx4)/6 + xt;
 yt = (ky1+2*ky2+2*ky3+ky4)/6 + yt;
 zt = (kz1+2*kz2+2*kz3+kz4)/6 + zt;

 xt_st(i) = xt;
 yt_st(i) = yt;
 zt_st(i) = zt;

 x_st(i) = x;
 y_st(i) = y;
 z_st(i) = z;

 x = x + h;
 y = y + h;
 z = z + h;

 if i < 51
     mprintf("x(t) = %f, y(t) = %f, z(t) = %f\n", xt_st(i), yt_st(i), zt_st(i));
 end;

end;

plot(x_st, xt_st, 'g*-');
plot(y_st, yt_st, 'b*-');
plot(z_st, zt_st, 'r*-');
h1 = legend(['X(t)'; 'Y(t)'; 'Z(t)']);
ylabel("t");

figure();
param3d(xt_st, yt_st, zt_st);
h1=legend(['X(t) = f(Y(t), Z(t))']);

//graphs after compression
minValue = min(xt_st);
maxValue = max(yt_st);
for i = 1:2000
 xt_st(i) = ((xt_st(i)-minValue) * 2)/(maxValue-minValue) - 1;
 yt_st(i) = ((yt_st(i)-minValue) * 2)/(maxValue-minValue) - 1;
 zt_st(i) = ((zt_st(i)-minValue) * 2)/(maxValue-minValue) - 1;
end
figure();

plot(x_st, xt_st,'g*-');
plot(y_st, yt_st,'b*-');
plot(z_st, zt_st,'r*-');
h1=legend(['X(t)'; 'Y(t)'; 'Z(t)']);
ylabel("t");

figure();
param3d(xt_st, yt_st, zt_st);