Textbook problem 6.13 Octave Script

% This file is the solution for Problem 6.13 in the text on induction machines.
%  See the notes from March 2 and following.  It demonstrates using inline
% fuctions in Gnu Octave, and using fsolve.
%
%  Rob Frohne and EMEC Class 2011
%
%
s=.035
R1=.21
X1=.26
X2=.26
Xm=10.1
V1=230/sqrt(3)
P=15000
nph=3
np=4
%x(1) = real(I2), x(2) = imag(I2), x(3) = Rr
a = @(x) [x(3)./s+x(3).*(1+X2./Xm)+X1.*x(3)./s./Xm];
b = @(x) [X1.*(1+X2./Xm)-R1.*x(3)./s./Xm+X2];
fcn = @(x) [x(1).*a(x)-x(2).*b(x)-V1, x(2).*a(x)+x(1).*b(x),x(1)-sqrt(181.5/x(3))];
[x,fval,info]=fsolve(fcn,[37,0,.13],optimset("ComplexEqn", "Off"));
I2=complex(x(1),x(2))
Rr=x(3)

wr=2*pi.*[-60:.01:800].*2./np;
wa=2*pi*60*2/np;
s=(wa-wr)./wa;

Il=V1./(R1+j*X1+(j*Xm.*(Rr./s+j*X2))./(Rr./s+j*X2+j*Xm));
I2=Il*j*Xm./(j*Xm+Rr./s+j*X2);
T=abs(I2.*I2).*Rr.*(1-s)./s./wr;
figure(1)
subplot(2,1,1)
plot(wr*60/2/pi,T)
title('Torque')
xlabel('rotor speed (RPM)')
subplot(2,1,2)
plot(s,T)
xlabel('slip, s')