pushover cureve by MATLAB

[hamkar]

Hi
is there any one to control my following MATLAB code to
pushover curve:

clear all;
clc;
d=45;
i=1;
%----------------------Distributed Plasticity Analysis---------------------
for w=0:20:300
ft=fopen('bstory.tcl','wt');
%units: kips.ft
fprintf(ft,'#--------------------------------------------------------#\n');
% fprintf(ft,'set dataDir saveRecorder \n\n');
fprintf(ft,'wipe \n');
fprintf(ft,'model basic -ndm 2 -ndf 3 \n');
fprintf(ft,'source Wsection.tcl \n');
fprintf(ft,'source LibUnits.tcl \n');
fprintf(ft,'node 1 0 0\n');
fprintf(ft,'node 2 0 [expr $cm*300]\n');
fprintf(ft, 'node 3 [expr $cm*600] [expr $cm*300] \n');
fprintf(ft, 'node 4 [expr $cm*600] 0 \n');
fprintf(ft,'#--------------------------------------------------------#\n');
fprintf(ft,'fix 1 1 1 1 \n');
fprintf(ft,'fix 4 1 1 1 \n');
fprintf(ft,'#--------------------------------------------------------#\n');
fprintf(ft, 'set Fy [expr 2400/101.973]; # Fy=2400 Kg/cm^2 \n');
fprintf(ft, 'set Es [expr 2100000/101.973]; # Steel Youngs Modulus \n');
fprintf(ft, '# ELEMENT properties \n');
fprintf(ft, '#uniaxialMaterial Steel02 $matTag $Fy $E $b $R0 $cR1 $cR2 $a1 $a2 $a3 $a4 \n');
fprintf(ft, 'uniaxialMaterial Steel02 1 33.1 289350 0.03 10 0.925 0.15 \n');
fprintf(ft,'geomTransf PDelta 1 \n');
fprintf(ft,'#--------------------------------------------------------#\n');
fprintf(ft,'# beam-Column sections: W14x53 \n');
fprintf(ft,'#set BeamCOLSecTag 1 \n');
fprintf(ft, 'set d [expr $cm*25]; # depth \n');
fprintf(ft, 'set bf [expr $cm*26]; # flange width \n');
fprintf(ft, 'set tf [expr $cm*1.25]; # flange thickness \n');
fprintf(ft, 'set tw [expr $cm*0.75]; # web thickness \n');
fprintf(ft, 'set nfdw 10; # number of fibers along dw \n');
fprintf(ft, 'set nftw 4; # number of fibers along tw \n');
fprintf(ft, 'set nfbf 10; # number of fibers along bf \n');
fprintf(ft, 'set nftf 4; # number of fibers along tf \n');
fprintf(ft, 'Wsection 1 1 $d $bf $tf $tw $nfdw $nftw $nfbf $nftf \n');
fprintf(ft,'#--------------------------------------------------------#\n');
fprintf(ft,'#element $eleTag $iNode $jNode $numIntgrPts $secTag \n');
fprintf(ft,'element nonlinearBeamColumn 1 1 2 5 1 1 \n');
fprintf(ft,'element nonlinearBeamColumn 2 2 3 5 1 1 \n');
fprintf(ft,'element nonlinearBeamColumn 3 3 4 5 1 1 \n');
fprintf(ft,'#--------------------------------------------------------#\n');
fprintf(ft,'recorder Node -file disp.out -time -node 2 -dof 1 disp \n');
fprintf(ft,'#--------------------------------------------------------#\n');
fprintf(ft,'set w %d \n',w);
fprintf(ft,'set b [expr 1.6*$w] \n');
fprintf(ft,'set c [expr 1.2*$d] \n');
fprintf(ft,'#--------------------------------------------------------#\n');
fprintf(ft,'pattern Plain 1 Constant { \n');
fprintf(ft, 'eleLoad -ele 2 -type -beamUniform -1; \n');
fprintf(ft,'} \n');
fprintf(ft,'pattern Plain 2 Linear { \n');
fprintf(ft,'load 2 $b 00.0 0.0 \n');
fprintf(ft,'} \n');
fprintf(ft,'#--------------------------------------------------------#\n');
fprintf(ft,'constraints Transformation \n');
fprintf(ft,'numberer RCM \n');
fprintf(ft,'system BandGeneral \n');
fprintf(ft,'test NormDispIncr 1.0e-6 6 \n');
fprintf(ft,'algorithm Newton \n');
fprintf(ft,'integrator LoadControl 0.01 \n');
fprintf(ft,'analysis Static \n');
fprintf(ft,'analyze 100 \n');
fprintf(ft,'loadConst -time 0.0 \n');
fprintf(ft,'#--------------------------------------------------------#\n');
fprintf(ft,'getTime \n');
fclose(ft);
%
!opensees bstory.tcl
D=zeros(100,1);
[step(:,1),D(:,i)]=textread('disp.out','%f %f');
maxD=max(abs(D(:,i)));
W(i,1)=w;
maxDD(i,1)=maxD;
i=i+1;
end
%
%---------------------------Plot w-deflection------------------------------
plot(maxDD,W);