function [t,y] = nareul(funname,timerange,x0)
% Narrated tutorial Euler method using same format as ode45
% funame    is a string containing the name of a function that gives the time
%           derivaative of the system you want to solve.  For more 
%           information on the format of this function, type 'help odefile'
% timerange is a two-element vector with start and end times
% x0        is the initial state of the system.
% y         gives the state of the system at times in vector
% t


% Charlie Sullivan, 6/00
% Modified 6/29/03:  minor changes to make it compatible with MATLAB 6.5
% (R13).
% Also fixed double-digit initial condition bug.

disp(' ')
disp('Welccome to nareul.  This program solves differential equations') 
disp('numerically using the basic, not terribly accurate Euler method.')
disp('It is intended as a tutorial, to show you how this works.  The highly')
disp('accurate MATLAB solvers such as ode45 use the same format as this')
disp('function, so once you learn to use this, you can use ode45 the same way')
disp('to get a more accurate, faster solution.')
disp('------------------------------------------------------------------------')

%string used to get linefeeds on screen
feed = [' ';' ';' ';' ';' ';' ';' ';' ';' ';' '];

% Check inputs:

if (nargin < 3)
disp('You have not included all the necessary inputs.')
disp('Please try again:  Call nareul with:')
disp('nareul(funname,timerange,x0)')
disp('where funnname is a string with the name of your derivative function,')
disp('t is a two-element vector with starting and ending times,')
disp('and x0 is a *column* vector with the initial state(s)')
return
end

if (nargin > 3)
disp('You have not included too many inputs.')
disp('nareul does not take any of the optional inputs that') 
disp('can be used with ode45')
disp('Please try again:  Call nareul with:')
disp('nareul(funname,timerange,x0)')
disp('where funnname is a string with the name of your derivative function,')
disp('t is a two-element vector with starting and ending times,')
disp('and x0 is a *column* vector with the initial state(s)')
return
end

if not(ischar(funname))
disp('Error:  Your function name was not a string')
disp('You can either set a variable equal to a string, and call')
disp('nareul using that string, e.g., to call derivfun.m')
disp('     funname = ''derivfun''')
disp('     nareul(funname,timerange,x0)')
disp('or you can put the string right in like:')
disp('     nareul(''derivfun'',t,x0)')
return
end

last_two_char = funname( (length(funname)-1):length(funname));
if (last_two_char=='.m')
disp('Error:  You do not need to include the .m when you tell nareul your')
disp('function name.')
return
end

if max(size(timerange) ~= [1 2])
disp('Error:  Your second argument must be a row vector containing')
disp('the initial and final times for which you want a solution')
disp('Example:  for times 0 to 12, you can use:')
disp('     times = [0 12])')
disp('     nareul(funname,times,x0)')
disp('or just')
disp('     nareul(funname,[0 12],x0)')
return
end

% Check that x0 is a column vector
x0size = size(x0);
if x0size(2) > 1
disp('Error:  Your third argument must be a column vector containing the initial')
disp('state.  (or a scalar if you have only one state.)')
disp('To create a column vector, use a semicolon between elements')
disp('Example:  for initial states x1(t=0) = 17 and x2(t=0) = -17:')
disp('     x0 = [0 ; -17])')
disp('     nareul(funname,timerange,x0)')
disp('or just')
disp('     nareul(funname,timrange,[0 ; -17])')
return 
end

% Initial Check on inputs is OK; now confirm with user:
disp('You have requested to solve the differential equations defined by the')
disp(['derivative function ', funname,'.m , for the time period from ' num2str(timerange(1))])
disp(['to ' num2str(timerange(2)) ', with the initial state'])
x0

contin = input('Hit enter to continue, 0 to abort   >');
if (length(contin) > 0) 
	disp('Aborted by user.  Start again.')
	return
end

%%%%%%%old code to use dialog box
%pushed = questdlg('OK?');

%if (length(pushed) ~= length('Yes')) 
%	disp('Aborted by user.  Start again.')
%	return
%end


%%%%%%% set up variables  %%%%%%%%%%%
n = length(x0);
t = timerange(1);
x = x0;
index =1; % number of point we are doing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Set up to plot
figure(1)
hold off
plot(timerange,zeros(size(timerange)),'.')
hold on
plot(timerange(1),x0,'o')
hold off
xlabel('time')
ylabel('states')
title('Initial conditions')

% Make a string that, executed, calls funnmae with proper parameters
x0str = num2str(x0);
string = ['xdot = ', funname, '(',num2str(timerange(1)),', ['];
for i = 1:(n-1)
	string = [string, x0str(i,:), '; '];
end
string = [string, x0str(n,:), ']);'];

disp(feed)
disp('Your initial conditions are now displayed in Figure 1')
disp(['Your derivative function, ',funname,'.m will now be called'])
disp('to calculate the derivative at the initial conditions,')
disp('using the syntax:')
disp(' ')
disp(string)
disp(' ')
disp('where the first argument is the time (in this case the initial time)')
disp('and the second argument is the state vector (in this case the state vector')
disp('at the initial conditions).')
disp('If the program crashes now, it is probably an error in your')
disp(['derivative function, ',funname,'.m.  You can troubleshoot'])
disp('by calling your function with the line above and making sure')
disp('it works.')
disp(' ')

contin = input('Hit enter to continue, 0 to abort   >');
if (length(contin) > 0) 
	disp(' ')
	disp(' ')
	disp('Aborted by user.  Start again.')
	return
end

%do it--actually calculate derivative
eval(string);

% Check xdot

xdotsize = size(xdot);

% Check that xdot is a column vector
if xdotsize(2) > 1
disp(feed)
disp('Error:  Your function must return the derivative in the form of a column vector')
disp('when you call it with the line')
disp(' ')
disp(string)
disp(' ')
disp('where the first argument is the time (in this case the initial time)')
disp('and the second argument is the state vector (in this case the state vector')
disp('at the initial conditions).')
disp('Please fix your function and try nareul again after your function is fixed.')
return 
end

% Check that xdot is the right length
if xdotsize(1) ~= n
disp(feed)
disp('Error:  Your function returned a vector of the wrong size.  It must') 
disp('return the derivative in the form of a column vector, the same size')
disp('as your initial condition vector:')
x0
disp('when you call your function with the line')
disp(' ')
disp(string)
disp(' ')
disp('where the first argument is the time (in this case the initial time)')
disp('and the second argument is the state vector (in this case the state vector')
disp('at the initial conditions).')
disp('Please fix your function and try nareul again after your function is fixed.')
return 
end

% %%%   Plot derivative   %%%
% tlength = (timerange(2)-timerange(1))/10;
% xvector = tlength * xdot;
% 
% 
% hold on
% quiver(t*ones(n,1),x(:,index),tlength*ones(n,1),xvector,0)
% hold off
% title('Derivatives at initial condition calculated by your function')

tlength = (timerange(2)-timerange(1))/20;

%%%   Plot derivative   %%%
%tlength = (timerange(2)-timerange(1))/10;
xvector = tlength * xdot;
		
hold on
quiver(t(index)*ones(n,1),x(:,index),tlength*ones(n,1),xvector,0)
hold off
title('Derivatives at initial condition calculated by your function')

% Enter loop to do more points

while max(t) < timerange(2)
	

	%msgbox('Derivatives calculated are shown.   Click on a time to calculate the new state at that time')
	disp(feed)
	disp('Derivatives calculated are shown.   Click on a time to calculate the')
	disp(['new state at that time.'])
	[tnew,xjunk] = ginput(1);
	if isempty(tnew), 
		tnew = t(index)+tlength;
		disp(feed)
		disp('No time selected; using default point spacing')
		disp('Next time, click on a time on the plot window to select manually')
		disp('Hit any key to continue.')
		pause
	end
	%Actually calculate new point now!
	dt = (tnew-t(index));
	x = [x (dt*xdot + x(:,index))];
	t = [t t(index)+dt];
	index = index + 1;
	
	%keyboard
	plot(timerange,zeros(size(timerange)),'.')
	hold on
	plot(t,x','o')
	plot(t,x')
	%t
	%x
	%dt
	%xdot
	hold off
	xlabel('time')
	ylabel('states')
	title('New points calculated')
	
	disp(feed)
	disp('The new point at the time you clicked has been calculated.')
        if max(t) > timerange(2)
              disp('Solution for requested time range completed.')
              disp('Nareul will now quit.')
              break
        end
	disp(['Your derivative function, ',funname,'.m will now be called'])
	disp(['again to calculate the derivative at this new time and state.'])
	contin = input('Hit enter to continue, 0 to abort   >');
	if (length(contin) > 0)
		disp(' ')
		disp(' ') 
		disp('Aborted by user.  Start again.')
		return
	end

	xdot = feval(funname,t(index), x(:,index));

	%%%   Plot derivative   %%%
	%tlength = (timerange(2)-timerange(1))/10;
	xvector = tlength * xdot;
	
	
	hold on
	quiver(t(index)*ones(n,1),x(:,index),tlength*ones(n,1),xvector,0)
	hold off
	title('Derivatives calculated by your function')
	
end

%return what we've been calling x as what ode45 calls y
y = x;

return
end

function big
return
end