clear;clc;
%output = fopen('output','w');
%global output2;
%output2 = fopen('output2','w');
output3 = fopen('output3','w');
%fprintf(output,'PART C\n');
%fprintf(output2,'PART D\n');
fprintf(output3,'PART E\n');
%fprintf(output2,'\n\t q \t h \t\t\t x \t\t\t max norm \t\t max norm h*|Phi-Phi_star|\n');
fprintf(output3,'\n\t q \t h \t\t\t x \t\t\t max norm \t\t max norm h^3*v\n');
for q = 0:1:3
    %fprintf(output,'\n\t q \t h \t\t\t max norm \t max |c2+s2-1|\n');
    
    figure(1);
    subplot(2,2,q+1)
    title(sprintf('q = %i',q))
    hold on 
    five_terms = 0;
    Tspan = [2*sqrt(q) 2*sqrt(q)+1];
    IV = [0 2*pi*sqrt(q) 1 0];
    h = 0.2;
    tolerance = 0.5e-6;
    [t,y,v] = Sys4ODEsRK4(@odes,Tspan,IV,h,q);
    %fprintf(output,'\t %i \t %f \t %d \t %f \n',q,h,max(...
        %[max(abs(y(3,:)-cos(pi/2*t.^2))) max(abs(y(4,:)-sin(pi/2*t.^2)))])...
        %,max(y(:,3).^2+y(:,4).^2-1));
    %while(max(abs(y(3,:)-cos(pi/2*t.^2))) && max(abs(y(4,:)-sin(pi/2*t.^2))) >= tolerance);
    for h = [.2 .025 .0125]
        h = h/2;
        [t,y,v] = Sys4ODEsRK4(@odes,Tspan,IV,h,q);
%        fprintf(output,'\t %i \t %f \t %d \t %f \n',q,h,max(...
%        [max(abs(y(3,:)-cos(pi/2*t.^2))) max(abs(y(4,:)-sin(pi/2*t.^2)))])...
%        ,max(y(:,3).^2+y(:,4).^2-1));
    end
     
    plot(t,y(3,:),'blue');plot(t,y(4,:),'Color',[0 0.5 0]);
    %plot(t,y(1,:),'cyan'):plot(t,y(2,:),'magenta');
        x = linspace(2*sqrt(q),2*sqrt(q)+1,200);
        c = cos(pi/2*x.^2);
        s = sin(pi/2*x.^2);
    %d2cdx2 = -pi^2*x.^2.*c-pi*s./sqrt(c.^2+s.^2);
    %d2sdx2 = -pi^2*x.^2.*c+pi*c./sqrt(c.^2+s.^2);
    plot(x,c,'--c');plot(x,s,'--g');
    %plot(t,d2cdx2,'--c');plot(t,d2sdx2,'--m');
%    exact = [c;s;dudx;dvdx];
    legend('u','v','c exact','s exact')
    hold off
     fig = figure(1);
     print(fig,'-dpsc','partd')
end
