function main()
   disp("Hello World");
end%This is a test case for solve initial value problem using euler's method

%function float output=dy(float t, float y)
%	%differential equation dy=(t-y)/2
%	output=(t-y)/2;
%end

%function float output=y(float t)
%	%Exact solution of dy
%	output=3*exp(-t/2)-2+t;
%end

function float output=dy(float t, float y)
	%differential equation 
	output=t^2-y;
end



function float output=y(float t)
	%Exact solution of d
	output=-exp(-t)+t^2-2*t+2;
end
function matrix E=Euler(float a, float b, float ya, int N)
	%Solve differential equation dy using Euler's method
	%y_k+1=y_k+hf(t_k, y_k) for k=0,1,2..N-1
	%a and b are the left and right end points
	%ya is the initial condition y(a)
	%N is the numbe rof steps
	%Output E, first column is abscissas and second coordinate is ordinates

	float h;


	int i;
	E=newmat(N+1,2);
	
	h=(b-a)/N;
	E[0,0]=a;
	for i=1:N
		E[i,0]=E[i-1,0]+h;
	end

	E[0,1]=ya;
	for i=0:N-1
		E[i+1,1]=E[i,1]+h*dy(E[i,0],E[i,1]);
	end

end


function main()
	%this function solves the differential equation dy=(t-y)/2 using euler's method and displays the Euler solution with the exact solution and percentage error

	matrix Reuler; %Solution using euler method
	matrix Moutput; %output matrix
	float Rexact; %exact solution
	int n;
	int i;

	Reuler=Euler(0., 2., 1., 50);

	n=height(Reuler);
	disp(n);
	Moutput=newmat(n,4);
	for i=0:n-1
		Rexact=y(Reuler[i,0]);
		Moutput[i,0]=Reuler[i,0]; %t
		Moutput[i,1]=Reuler[i,1]; %Euler method solution
		Moutput[i,2]=Rexact; %exact solution
		Moutput[i,3]=(Moutput[i,2]-Moutput[i,1])/Moutput[i,2]*100.0; % percent error
	end
	disp("        t         Euler       Exact       %Error");
	disp(Moutput);
	save(Moutput,"MoutputEuler");

end%This is a test case for solve initial value problem using Heun's method

function float output=dy(float t, float y)
	%differential equation dy=(t-y)/2
	output=(t-y)/2;
end

function float output=y(float t)
	%Exact solution of dy
	output=3*exp(-t/2)-2+t;
end
function matrix E=Heun(float a, float b, float ya, int N)
	%Solve differential equation dy using Heun's method
	%a and b are the left and right end points of the independent variable
	%ya is the initial condition y(a) of the dependent variable
	%N is the numbe rof steps
	%Output E, first column is indepnedent variable and second variable is dependent variable

	float h;
	float k1;
	float k2;

	int i;
	E=newmat(N+1,2);
	
	h=(b-a)/N;
	E[0,0]=a;
	for i=1:N
		E[i,0]=E[i-1,0]+h;
	end

	E[0,1]=ya;
	for i=0:N-1
		k1=dy(E[i,0],E[i,1]);
		k2=dy(E[i+1,0], E[i,1]+h*k1);
		E[i+1,1]=E[i,1]+h/2*(k1+k2);
	end

end


function main()
	%this function solves the differential equation dy=(t-y)/2 using Heun's method and displays the solution with the exact solution and percentage error

	matrix Rheun; %Solution using Heun's method
	matrix Moutput; %output matrix
	float Rexact; %exact solution
	int n;
	int i;

	Rheun=Heun(0., 2., 1., 50);

	n=height(Rheun);
	disp(n);
	Moutput=newmat(n,4);
	for i=0:n-1
		Rexact=y(Rheun[i,0]);
		Moutput[i,0]=Rheun[i,0]; %t
		Moutput[i,1]=Rheun[i,1]; %Euler method solution
		Moutput[i,2]=Rexact; %exact solution
		Moutput[i,3]=(Moutput[i,2]-Moutput[i,1])/Moutput[i,2]*100.0; % percent error
	end
	disp("        t         Euler       Exact       %Error");
	disp(Moutput);
	save(Moutput,"MoutputHeun");

end%This is a test case for solving initial value problem using 4th order Runge-Kutta method

%function float output=dy(float t, float y)
%	%differential equation dy=(t-y)/2
%	output=(t-y)/2;
%end

function float output=dy(float t, float y)
	%differential equation 
	output=t^2-y;
end

%function float output=y(float t)
%	%Exact solution of dy
%	output=3*exp(-t/2)-2+t;
%end

function float output=y(float t)
	%Exact solution of dy
	output=-exp(-t)+t^2-2*t+2;
end

function matrix E=RK4(float a, float b, float ya, int N)
%solves the diff eq using 4th order runge kutta method
%a is the initial value of the independent variable
%b is the end value of the indepenedent variable
%ya is the initial value of the the dependent variable
%N is the number of steps
%output is a matrix with independent variable in the first column and dependent variable in the second column

	float h;
	float k1;
	float k2;
	float k3;
	float k4;
	int i;
	
	E=newmat(N+1,2);
	h=(b-a)/N;
	E[0,0]=a;
	
	for i=1:N
		E[i,0]=E[i-1,0]+h;
	end

	E[0,1]=ya;
	for i=0:N-1
		k1=h*dy(E[i,0],E[i,1]);
		k2=h*dy(E[i,0]+h/2., E[i,1]+k1/2.);
		k3=h*dy(E[i,0]+h/2., E[i,1]+k2/2.);
		k4=h*dy(E[i,0]+h,E[i,1]+k3);
		E[i+1,1]=E[i,1]+(k1+2*k2+2*k3+k4)/6.0;
	end

end


function main()
	%this function solves the differential equation dy using 4th order Runge Kutta method and displays the RK4 solution with the exact solution and percentage error

	matrix R; %Solution using RK4 method
	matrix Moutput; %output matrix
	float Rexact; %exact solution
	int n;
	int i;

	
	R=RK4(0., 2., 1., 50);

	n=height(R);

	Moutput=newmat(n,4);
	for i=0:n-1
		Rexact=y(R[i,0]);
		Moutput[i,0]=R[i,0]; %t
		Moutput[i,1]=R[i,1]; %Euler method solution
		Moutput[i,2]=Rexact; %exact solution
		Moutput[i,3]=(Moutput[i,2]-Moutput[i,1])/Moutput[i,2]*100.0; % percent error
	end
	disp("        t              RK4        Exact       %Error");
	disp(Moutput);
	save(Moutput,"MoutputRK4");

end%This is a test case for solving 2nd order differental equaion initial value problem using 4th order Runge-Kutta method

%The diff eq is: x''(t)+4x'(t)+5x(t)=0
%reformulate to the following form for solution:
%dx/dt=y, dy/dt=-5x-4y
%initial conditions: x(0)=3, y(0)=-5

%function float output=dy(float x, float y)
%	%differential equation 
%	output=-5*x-4*y;
%end
function float output=dy(float x, float y)
	output=3*x+2*y;
end

function float output=dx(float x, float y)
	output=x+2*y;
end

%function float output=x(float t)
%	%Exact solution of dx
%	output=3*exp(-2*t)*cos(t)+exp(-2*t)*sin(t);
%end
function float output=x(float t)
	%Exact solution of dx
	output=4*exp(4*t)+2*exp(-t);
end

function matrix E=RK4_2(float a, float b, float ya, float xa, int N)
%solves the diff eq using 4th order runge kutta method
%a is the initial value of the independent variable
%b is the end value of the indepenedent variable
%ya is the initial value of dy/dt
%xa is the initial value of dx/dt

%N is the number of steps
%output is a matrix with independent variable in the first column and dependent variable in the second column

	float h;
	float g1;
	float g2;
	float g3;
	float g4;
	float f1;
	float f2;
	float f3;
	float f4;
	float y0;
	float y;
	int i;
	
	E=newmat(N+1,2);
	h=(b-a)/N;
	E[0,0]=a;
	E[0,1]=xa;
	
	
	for i=1:N
		E[i,0]=E[i-1,0]+h;
	end

	y0=ya;
	disp(y);
	for i=0:N-1
		g1=dy(E[i,1],y0);
		f1=dx(E[i,1],y0);
		
		g2=dy(E[i,1]+h/2*f1, y0+h*g1/2.);
		f2=dx(E[i,1]+h/2*f1, y0+h/2.0*g1);
		
		g3=dy(E[i,1]+h/2*f2, y0+h*g2/2.);
		f3=dx(E[i,1]+h/2*f2, y0+h/2.0*g2);
		
		g4=dy(E[i,1]+h*f3,y0+h*g3);
		f4=dx(E[i,1]+h*f3, y0+h*g3);
		
		y=y0+h*(g1+2*g2+2*g3+g4)/6.0;
		
		
	
		
		
		E[i+1,1]=E[i,1]+h*(f1+2*f2+2*f3+f4)/6.0;
		y0=y;
	end

end


function main()
	%this function solves the differential equation dy using 4th order Runge Kutta method and displays the RK4 solution with the exact solution and percentage error

	matrix R; %Solution using RK4 method
	matrix Moutput; %output matrix
	float Rexact; %exact solution
	int n;
	int i;
	setFormat("long");
	R=RK4_2(0.0, 0.2, 4.0, 6.0, 10);

	n=height(R);
	disp(n);
	Moutput=newmat(n,4);
	for i=0:n-1
		Rexact=x(R[i,0]);
		Moutput[i,0]=R[i,0]; %t
		Moutput[i,1]=R[i,1]; %Euler method solution
		Moutput[i,2]=Rexact; %exact solution
		Moutput[i,3]=(Moutput[i,2]-Moutput[i,1])/Moutput[i,2]*100.0; % percent error
	end
	disp("        t         Euler       Exact       %Error");
	disp(Moutput);
	save(Moutput,"MoutputRK4_2");

end%This is a test case for solving 2nd order differental equaion initial value problem using 4th order Runge-Kutta method

%The diff eq is: x''(t)+4x'(t)+5x(t)=0
%reformulate to the following form for solution:
%dx/dt=y, dy/dt=-5x-4y
%initial conditions: x(0)=3, y(0)=-5

function float output=dy(float x, float y)
	%differential equation 
	output=-5*x-4*y;
end

function float output=dx(float x, float y)
	output=y;
end

function float output=x(float t)
	%Exact solution of dx
	output=3*exp(-2*t)*cos(t)+exp(-2*t)*sin(t);
end

%function float output=dy(float x, float y)
%	output=3*x+2*y;
%end

%function float output=dx(float x, float y)
%	output=x+2*y;
%end

%function float output=x(float t)
%	%Exact solution of dx
%	output=4*exp(4*t)+2*exp(-t);
%end

function matrix E=RK4_2(float a, float b, float ya, float xa, int N)
%solves the diff eq using 4th order runge kutta method
%a is the initial value of the independent variable
%b is the end value of the indepenedent variable
%ya is the initial value of dy/dt
%xa is the initial value of dx/dt

%N is the number of steps
%output is a matrix with independent variable in the first column and dependent variable in the second column

	float h;
	float g1;
	float g2;
	float g3;
	float g4;
	float f1;
	float f2;
	float f3;
	float f4;
	float y0;
	float y;
	int i;
	
	E=newmat(N+1,2);
	h=(b-a)/N;
	E[0,0]=a;
	E[0,1]=xa;
	
	
	for i=1:N
		E[i,0]=E[i-1,0]+h;
	end

	y0=ya;
	disp(y);
	for i=0:N-1
		g1=dy(E[i,1],y0);
		f1=dx(E[i,1],y0);
		
		g2=dy(E[i,1]+h/2*f1, y0+h*g1/2.);
		f2=dx(E[i,1]+h/2*f1, y0+h/2.0*g1);
		
		g3=dy(E[i,1]+h/2*f2, y0+h*g2/2.);
		f3=dx(E[i,1]+h/2*f2, y0+h/2.0*g2);
		
		g4=dy(E[i,1]+h*f3,y0+h*g3);
		f4=dx(E[i,1]+h*f3, y0+h*g3);
		
		y=y0+h*(g1+2*g2+2*g3+g4)/6.0;
		
		
		E[i+1,1]=E[i,1]+h*(f1+2*f2+2*f3+f4)/6.0;
		y0=y;
	end

end


function main()
	%this function solves the differential equation dy using 4th order Runge Kutta method and displays the RK4 solution with the exact solution and percentage error

	matrix R; %Solution using RK4 method
	matrix Moutput; %output matrix
	float Rexact; %exact solution
	int n;
	int i;
	
	R=RK4_2(0.0, 5.0, -5.0, 3.0, 50);
	setFormat("short");
	n=height(R);
	disp(n);
	Moutput=newmat(n,4);
	for i=0:n-1
		Rexact=x(R[i,0]);
		Moutput[i,0]=R[i,0]; %t
		Moutput[i,1]=R[i,1]; %Euler method solution
		Moutput[i,2]=Rexact; %exact solution
		Moutput[i,3]=(Moutput[i,2]-Moutput[i,1])/Moutput[i,2]*100.0; % percent error
	end
	disp("        t              RK4        Exact       %Error");
	disp(Moutput);
	save(Moutput,"MoutputRK4_2");

endfunction  main()
%test for loop
  int x;
  int i;
  i=0;
  for i=10:-1:0
      x=x+1;
  end
  disp(x);
end
function main()
%test nested loop
  float i;
  float j;
  matrix a;
  a=newmat(3,3);
  for i=1.0:3.0
    for j=1.0:3.0
      disp(i+j);
    end
  end

  
end
function  main()

%nested for increments

  int i;
  int j;

  for i=1:1:3
    for j=1:1:3
      disp(i*j);
    end
  end

  
end
function  main()

%nested for increments

  float i;
  float j;

  for i=0.0:0.2:0.4
    for j=0.0:0.2:0.4
      disp(i+j);
    end
  end

  
end
function  main()

%nested for increments

  float i;
  float j;

  for i=0.4:-0.2:0.0
    for j=0.4:-0.2:0.0
      disp(i-j);
    end
  end

  
end
function main()
%test while with for
  int i;
  int j;
  i = 5;
  
  while i > 0 
    
    i = i - 1;
    for j=0:3
        disp(i+j);
    end
  end
  
  disp("i="+i);
  
end
% OK! call a function with return type without assinging it to other variable  
function int m = test (int x, int y, int z)
    
    m=x+y+z;   

end

function  main ()
   
  
   int x;
   int y;
   int z;
   x=1;
   y=2;
   z=3;
   disp(test(x,y,z));
   
end
function main()
int x;
int y;
x=3;
y=4;
    if x==3
       x=x+1;
    else
       x=x-1;
    end   
    disp(x);
end
function main()
int x;
int y;
x=3;
y=4;
    if x==3
       if y==4 
          x=x+1;
       else 
          x=x-1;
       end
    else
       if x==5
          x=x-1;
       else
       end
    end 
    disp(x);
end

function  main()
int x;
int y;
x=3;
y=4;
    if x==3
       
    else
       if y==4 
          x=x+1;
       else 
          if y==4 
            x=x+1;
              if y==5
                x=x+1;
              else
              end
          else 
            x=x-1;
          end
       end
    end
    disp(x);      
end
function main()
int x;
int y;
x=3;
y=0;

    if x==4
      y=4;
    elseif x==5
      y=5;
    elseif x==6
      y=6;
    elseif x==3
      y=3;
    
    else
    end
      
disp(y);  
       
end
function  main()
int x;
int y;
x=3;
y=0;

    if x==4
      y=4;
    elseif x==5
      y=5;
    elseif x==6
      y=6;
    elseif x==3
       if y==0
          y=y+0;
       else
          y=y+3;
       end   
    else
    end
      
disp(y);  
       
end
function main()
int x;
int y;
x=3;
y=0;

    if x==4
      y=4;
    elseif x==5
      y=5;
    elseif x==6
      y=6;
    elseif x==3
      y=3;
    else
        if x==4
           y=4;
        elseif x==5
           y=5;
         else
        end
    
    end
      
disp(y);  
       
end

function main()
%test simple addition
   float r;
   r=3.2+2.2;
   disp("r=3.2+2.2="+tostring(r));
endfunction main()
%test simple addition with negative numbers
   float r;
   r=-3.2+-2.2;
   disp("r=3.2+2.2="+tostring(r));
endfunction main()
%test floating point math with e
   float r;
   r=-3.+10.+1e3;
   disp("r=-3+10+1e3="+tostring(r));
endfunction main()
%test a complex chain of operations with parenthesis
   float r;
   r=(1.0+2.)*3.0;
   r=2.^2.;
   r=1.+1.2*1.2-2/2*4.0+3^2^2+4.^(1./2.)+sqrt(9.)*(1.+2.-3./6.);
   setFormat("long");
   disp("r=");
   disp(r);
endfunction main()
%test a complex chain of operations with parenthesis, log, exponential and trig functions
   float r;
   float pi;
   pi=3.1415926536;
   r=abs(-1.0)*log(exp(2.0))*log10(10.)-(sin(pi)+cos(pi)+tan(pi));
   setFormat("shortE");
   disp("r=");
   disp(r);
endfunction main()
%testa series of trig functions and mod function, and nested sqrt with trig function.
   float r;
   float a;
   float b;
   float c;
   float pi;
   pi=3.1415926536;
   a=3.;
   b=4.;
   c=5.;
   r=acos(a/b);
   disp(r);
   r=asin(a/b);
   disp(r);
   r=atan(a/b);
   disp(r);
   r=atan2(b,a);
   disp(r);
   r=acos(mod(c,a)/c);
   disp(r);
   r=sqrt(sin(pi/2.));
   disp(r);
endfunction main()
%test matrix addition
matrix m;
matrix x;
m=[1,2;3,4;5,6];
x=[1,1;2,1;3,4];
disp(m+1);
disp("");
disp(1+m);
disp("");
disp(2.0+m);
disp("");
disp(m+2.0);
disp("");
disp(m+x);
endfunction main()
%test matrix assignment and for loop
	matrix m;
	int i;
	m=newmat(1,6);
	for i=0:5
		m[0,i]=i;
	end
	
	for  i=0:5
		disp(m[0,i]);
	end

end

function main()
%test matrix save and load
	matrix m;
	matrix s;
	m=[1,0,0,1;0,1,0,1;0,0,1,1];
	disp("width of m=" + tostring(width(m)));
	disp("height of m=" + tostring(height(m)));
	s=m[0:2,0:2];
	disp(inv(s));
		
end

function main()
%test matrix scaler division
matrix m;
matrix x;
m=[2,2;2,2];
x=[2;2];
disp(m/1);
disp("");
disp(1/m);
disp("");
disp(2.0/m);
disp("");
disp(m/2.0);

endfunction main()
%test matrix element by element operations
matrix m;
matrix x;
m=[1,2,3,4];
x=[2,4,6,8];
disp(m./x);
disp("");
disp(m.\x);
disp("");
disp(m.*x);
disp("");
disp(m^2.0);
disp("");
disp(x^1);

endfunction main()
%test matrix addition
matrix m;
matrix x;
m=[1,1;1,1];
x=[2;2];
disp(m*1);
disp("");
disp(1*m);
disp("");
disp(2.0*m);
disp("");
disp(m*2.0);
disp("");
disp(m*x);
endfunction main()
%test matrix range access
	matrix m;
	matrix s;
	int i;
	m=newmat(2,6);
	for i=0:5
		m[0,i]=i;
	end

	for i=0:5
		m[1,i]=i*2;
	end
	
	s=m[0:1,3:4];
	disp("m=");
	disp(m);
	disp("m(0:1,3:4)=");
	disp(s);
	disp("m(:,:):");
	
end

float pi;

function main()
%test matrix operations by testing a rotate vector function
	matrix v;
	pi=3.1415926;
	v=[0,0,1]; %unit vector pointing in z direction
	disp(rotatex(v,30.0)); %rotate about x axis 30 deg, display resulting rotated vector
end

function matrix v = rotatex(matrix m, float angle)
	matrix Rx;
	Rx=[1,0,0;0,cos(angle*pi/180.0),-sin(angle*pi/180.0);0,sin(angle*pi/180.),cos(angle*pi/180.0)];
	v=m*Rx;
endfunction main()
%test matrix save and load
	matrix m;
	matrix s;
	int i;
	m=newmat(2,6);
	for i=0:5
		m[0,i]=i;
	end
	for i=0:5
		m[1,i]=i*2;
	end
	
	save(m,"msaveTest");
	s=load("msaveTest");
	disp(s);
		
end

function main()
%test matrix transpose
matrix m;
m=[1,2;3,4;5,6];
disp(m);
disp("");
disp(m');
endfunction  main()
   int i;
   int j;
   int k;
   i = 3;
   j = ---6;
   
   k = - 5 - - j + -3;
   disp(k);
   
end
%global variable name is the same with global function name
int x;
function int m=x()
   m=2;
end
function  main()
      disp(x);
      disp(x());
end
%local variable is the parameter of the global function
function int x=x()
   x=2;
end
function  main()
      disp(x());
end
%local variable is inside the scope of the global function
function int m=x()
  int x;
   x=2;
   m=x;
end
function  main()
      disp(x());
end

%local variable is inside the scope of the global function
function int m=test()
  int x;
   x=2;
   m=x;
end
function  main()
   int test;
   test=7;
   disp(test());
   disp(test);
end
int x;
function int x=test()
  x=4;
end
function main()
 
   disp(x);
   disp(test());
  
end
int x;
function int m=test()
  int x;
  x=3;
  m=x;
end
function  main()
 
   disp(x);
   disp(test());
  
end% OK -1 nested recursive
function int m = foo (int x)
    if(x > 0)
        m=foo(x-1);
    else 
        m=-1;
    end
end

function int m = bar (int x)
    m=x+1;
end

function main() 
    disp(foo(bar(5)));
end
% OK -1 nested recursive
function int m = foo (int x)
    if(x > 0)
        m=foo(x-1);
    else 
        m=-1;
    end
end

function int m = bar (int x)
    m=x+1;
end

function  main()
    disp(foo(bar(foo(5))));
end
