MP5.3_Runge-Kutta

//Machine Problem No. 5.3 - 4th-order Runge-Kutta Method
//code written by RCQuilala - MSCE (S)

#include <iostream>
#include <cmath>

using namespace std;

double fy(double t, double y, double z)
{
    double p=sin(t)+cos(y)+sin(z);
    return p;
}

double fz(double t, double z)
{
    double p=cos(t)+sin(z);
    return p;
}

int main()
{
    cout<<"Solve the system of ordinary differential equations"<<endl<<endl;
    cout<<"     dy/dt = sin t + cos y + sin z,   y(0) = 0"<<endl;
    cout<<"     dz/dt = cos t + sin z,           z(0) = 0"<<endl<<endl;
    cout<<"     t ϵ [0, 20] with 100 intervals,"<<endl<<endl;
    cout<<"using a Fourth-order Runge-Kutta Method."<<endl<<endl;
    cout<<" t           y           z"<<endl<<endl;

    double h=0.2, k1, k2, k3, k4, l1, l2, l3, l4;
    int i=0, n=100;
    double t[n+1], y[n+1], z[n+1];

    t[0]=0, y[0]=0, z[0]=0;
    cout<<t[0]<<"       "<<y[0]<<"      "<<z[0]<<endl;
    
    for(i=1;i<n+1;i++)
    {
        t[i]=(i)*h;
        
        k1=h*fy(t[i],y[i-1],z[i-1]);
        l1=h*fz(t[i],z[i-1]);
        k2=h*fy(t[i]+h/2,y[i-1]+k1/2,z[i-1]+l1/2);
        l2=h*fz(t[i]+h/2,z[i-1]+l1/2);
        k3=h*fy(t[i]+h/2,y[i-1]+k2/2,z[i-1]+l2/2);
        l3=h*fz(t[i]+h/2,z[i-1]+l2/2);
        k4=h*fy(t[i]+h,y[i-1]+k3,z[i-1]+l3);
        l4=h*fz(t[i]+h,z[i-1]+l3);
        
        y[i]=y[i-1]+(k1+2*k2+2*k3+k4)/6;
        z[i]=z[i-1]+(l1+2*l2+2*l3+l4)/6;
        
        cout<<t[i]<<"      "<<y[i]<<"       "<<z[i]<<endl;
    }
    
    return 0;
}