function [soln,dsoln,ddsoln] = calcshi(n)
% programme to plot the hrr stress distribution wrt angular coordinate \theta
% Mode I case
% This is a programme to calculate shi(theta) as explained
% in the class and in the assignment sheet. This programme will
% not run unless you fill in the blanks marked by '******'. To do that
% you need to understand the structure of the programme and know what
% you want to do. Taking a look at the Matlab ode solvers, like ode45,
% ode23s will help. Just type 'help ode45' at the matlab prompt.
% Once the blanks are filled in, it can be run by typing 
% [shi dshi ddshi]=calcshi(n);, where n is the value of the work hardening
% exponent that you have to use.
global N
 
N=1/n;
tol=0.0001;
ddshiatzero=input('Second deriv at zero?');
indicator=0;
while indicator == 0
y0=*******;
[t,shi]=ode45('hrr',pi/100:pi/50:pi,y0);
% check values at pi
nn=size(shi,1);
shiatpi=******;
dshiatpi=******;
plot(t,******);
hold on;
if abs(shiatpi) <= tol
   indicator=1;
   disp('Convereged within the given tolerance');
   soln=******;
   dsoln=******;
   ddsoln=******;
else
   indicator=0;
   str1='Value of shi at pi =';
   str2='Value of first derivative of shi at pi =';
   str3=num2str(shiatpi);
   str4=num2str(dshiatpi);
   str5=[str1 str3];
   str6=[str2 str4];
   disp(str5);
   disp(str6);
   ddshiatzero=input('Second deriv at zero?');
   
end   
end
