%clear all;
%close all;
%                  Program mtc (method of characteristic)
%
%     This program solves the initial value problem 
%
%            u_t + uu_x = 0,  u(x,0) = f(x)
%
%   using the method of characteristics.  There are four choices of
%   initial data. Choices 1 and 2 should develop a shock, but instead
%   'roll over' into a curve which no longer defines a function.
%    Choice 3 is a rarefaction wave.  The initial data for choice 4
%   is a small perturbation of the constant function f(x) = 1. In this
%   case the program computes the solution of the nonlinear equation 
%   and the solution of the linearized problem which is u_t + u_x = 0.
%       At run time user enters the choice of data, 1,2,3, or 4, and
%   the times t1,t2,t3,and t4 at which the solution is to be viewed.
%   The solution profiles at times t = 0,t1,t2,t3,t4 are all plotted on
%   the interval [-1 6]. The solution profiles of both linear and nonlinear
%   problems are plotted together in the case of data option 4.


disp('  ')
disp('Enter 1 for decreasing profile  ')
disp('Enter 2 for another decreasing profile  ')
disp('Enter 3 for increasing profile  ')
disp('Enter 4 for the small perturbation of f(x)= 1 ')
disp('Enter 5 for po')
m = input( 'enter the choice of initial data: 1,2,3,4 or 5   ')

t = input(' Enter the times in the form  [t1 t2 t3 t4 ]  ')


x = -5:.05:6;
n = length(x);

if m == 1

   ymiddle = 1.0 - .125*x.^2.*(3.0 - x);

  y = (x < 0) + ( (x < 2.001 ) - (x< 0)  ).*ymiddle  + .5*(x > 2);


elseif m == 2

     ymiddle  = .25*(x+1).*(x-2).^2;

   y = (x < 0) + ( (x < 2.001 ) - (x < 0) ).*ymiddle ;

elseif m == 3

ymiddle =  .5 - (x+1).^2 .*(x-2)/8 ;

  y =  .5*(x < -1) + ( ( x < 1.001 ) - ( x < -1) ).*ymiddle + (x > 1);

elseif m == 4

   ymiddle = .05*x.^3 - .15*x.^2 + 1.1;

    y = 1.1*(x < 0) + ( (x < 2.001) - (x < 0) ) .*ymiddle + .9*(x > 2);
    
elseif m==5
    
    %y=pi/2+atan(3*x);
    %y=pi/2-atan(x);
    
    %ymiddle=1-0.5*x;
    %y = (x < 0)*1 + ( (x < 1.001 ) - (x< 0)  ).*ymiddle  + .5*(x > 1);
    
    ymiddle=1-.5*x;
    y = (x < 0)*1 + ( (x < 1.001 ) - (x< 0)  ).*ymiddle  + 0.5*(x > 1);

end

   high = max(y); low = min(y); gap =.25*(high - low);
   x1 = x + t(1)*y;
   x2 = x + t(2)*y;
   x3 = x + t(3)*y;
   x4 = x + t(4)*y;
   figure
   plot(x,y, 'r')
   axis([-5,6, low-gap, high+gap])
   hold on
   
   plot(x1,y,'b-o')   
   plot(x2,y,'y-+')
   plot(x3,y,'k-.')
   plot(x4,y,'g--') 

hold off

say1=['t=' num2str(t(1))];
say2=['t=' num2str(t(2))];
say3=['t=' num2str(t(3))];
say4=['t=' num2str(t(4))];

legend('t=0',say1,say2,say3,say4,'Location', 'Southeast' )
   
if m == 4

  xl1 = x + t(1);
  xl2 = x + t(2);
  xl3 = x + t(3);
  xl4 = x + t(4);


  plot(xl1,y,'g')
  plot(xl2,y,'g')
  plot(xl3,y,'g')
  plot(xl4,y,'g')
else

end

%set(gca, 'XLim',[-1 10])

hold off
%legend('1','2','3','4','5')

