//  Rungen2.java: 2元連立常微分方程式の4次Runge-Kutta法による解法

import java.io.*;

public class Rungen2{
 
    public static void main( String argv[] ) throws IOException {

	double min = 0.0;
	double max = 10.0;
    
	double x, y0, y1, dx, h, yreal, yerror, maxerr;
	double t10, t11, t20, t21, t30, t31,       // Runge-Kutta法のための 
	    k10, k11, k20, k21, k30, k31 ,k40, k41;  // 一時変数        

	String s;   

	System.out.println("Input: dx");
	BufferedReader buf = new BufferedReader(new InputStreamReader(System.in));
	s = buf.readLine();
	dx = Double.parseDouble(s);
	System.out.println("dx=" +dx);

	h=dx/2.0;
   
	x = min;
	y0 = 0.0;                                // yの初期条件
	y1 = 1.0;                                // dy/dxの初期条件
	yreal = y0;
	yerror = 0.0;
	maxerr=0.0;

	try {
	    PrintWriter pw = new PrintWriter(new FileWriter("rungen.dat"));
	    
	    pw.println(x+" "+y0+" "+yreal+" "+yerror);

	    while(x <= max)
		{
		    
		    k10 = dx*f(x, y0, y1, 0);
		    t10 = y0+0.5*k10;
		    k11 = dx*f(x, y0, y1, 1);
		    t11 = y1+0.5*k11;
		    
		    k20 = dx*f(x+h, t10, t11, 0);
		    t20 = y0+0.5*k20;
		    k21 = dx*f(x+h, t10, t11, 1);
		    t21 = y1+0.5*k21;
        
		    k30 = dx*f(x+h, t20, t21, 0);
		    t30 = y0+    k30;
		    k31 = dx*f(x+h, t20, t21, 1);
		    t31 = y1+    k31;
        
		    k40 = dx*f(x + dx, t30, t31, 0);
		    k41 = dx*f(x + dx, t30, t31, 1);
        
		    y0 += (k10+2*k20+2*k30+k40)/6.0;
		    y1 += (k11+2*k21+2*k31+k41)/6.0;
		    
		    x += dx;    
		    yreal=Math.sin(x);              // 解析解
		    yerror=y0-yreal;                // 誤差=数値解-解析解    

		    if(Math.abs(yerror) > maxerr){
			maxerr=Math.abs(yerror);    // 誤差の絶対値の最大値  
		    }

		    //		    System.out.println("x= " +x + " y= " +y0);

		    pw.println(x+" "+y0+" "+yreal+" "+yerror);

		}
      
	    pw.close();

	    System.out.println("maxerr= " +maxerr);
	    System.out.println("data stored in rungen.dat");
      
	} catch (Exception e) {
	    System.err.println("File output error: " + e);
	}

        
    }

    public static double f(double x, double y0, double y1, int i)
    {
	if (i == 0) return(y1);               // d(y0)/dx 
	if (i == 1) return(-y0);              // d(y1)/dx 
	return((double)0.0);
    }
    
}