/* ode.c -- Solve initial value problem using explicit methods
   Written Jan 1, 2008 by Eric Olson for Mathematics 467

   For more information see Step 33 in First Steps in Numerical
   Analysis, 2nd Edition, by Hosking, Joe, Joyce and Turner */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double f(double y,double t){
    return t*t*(y+1);
}

double euler(int n, double y0, double t0, double T){
    int i;
    double yn,h,tn;
    h=(T-t0)/n;
    yn=y0;
    for(i=0;i<n;i++){
        tn=t0+i*h;
        yn=yn+h*f(yn,tn);
    }
    return yn;
}

double RK2(int n, double y0, double t0, double T){
    int i;
    double yn,h,tn;
    double k1,k2;
    h=(T-t0)/n;
    yn=y0;
    for(i=0;i<n;i++){
        tn=t0+i*h;
        k1=h*f(yn,tn);
        k2=h*f(yn+k1,tn+h);
        yn=yn+(1.0/2.0)*(k1+k2);
    }
    return yn;
}

int main(){
    double y;
    printf("Solve initial value problem using explicit methods.\n");
    y=euler(100,.3,0.0,1.0);
    printf("f=t*t*(y+1)\n");
    printf("y(0.0)=0.3\n");
    printf("Euler y(1.0)=%g\n",y);
    y=RK2(100,.3,0.0,1.0);
    printf("RK2 y(1.0)=%g\n",y);
    return 0;
}