Getting a warning while running the program.Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.394922e-018..

I am running the below given program which gives the warning message.How to remove the warning. ---------------------------------------------------------------------

   A barrel vault has radius r=25 ft, subtended angle of 80 degrees,
   length 50 f, and thickness 3in.Elastic modulus is 3 msi,Poisson's ratio is 0,and weight of the shell is 90 lb per square ft. Two curved edges are
   supported by rigid diaphragm and the other two edges are free. 
   Solve using 4 by 4 elements of a quarter of the shell.         
 Variable descriptions                                                      
%   k = element matrix in the local axes                                            
%   ke = element matrix in the global axes
%   kb = element matrix for bending stiffness
%   ks = element matrix for shear stiffness
%   km = element matrix for membrane stiffness
%   f = element vector
%   kk = system matrix                                             
%   ff = system vector                                                 
%   disp = system nodal displacement vector
%   gcoord = coordinate values of each node
%   nodes = nodal connectivity of each element
%   index = a vector containing system dofs associated with each element     
%   pointb = matrix containing sampling points for bending term
%   weightb = matrix containing weighting coefficients for bending term
%   points = matrix containing sampling points for shear term
%   weights = matrix containing weighting coefficients for shear term
%   bcdof = a vector containing dofs associated with boundary conditions     
%   bcval = a vector containing boundary condition values associated with    
%           the dofs in 'bcdof'                                              
%   kinmtsb = matrix for kinematic equation for bending
%   matmtsb = matrix for material property for bending
%   kinmtss = matrix for kinematic equation for shear
%   matmtss = matrix for material property for shear
%   kinmtsm = matrix for kinematic equation for membrane
%   matmtsm = matrix for material property for membrane
%   tr3d = transformation matrix from local to global axes

%-------------------------- % input data for control parameters %---------------------------

 clear
 nel=16;                   % number of elements
 nnel=4;                  % number of nodes per element
 ndof=6;                  % number of dofs per node
 nnode=25;                 % total number of nodes in system
 sdof=nnode*ndof;         % total system dofs  
 edof=nnel*ndof;          % degrees of freedom per element
 emodule=3e6;            % elastic modulus
 poisson=0.0;             % Poisson's ratio
 t=3;                   % plate thickness
nglxb=2; nglyb=2;        % 2x2 Gauss-Legendre quadrature for bending 
nglb=nglxb*nglyb;        % number of sampling points per element for bending
nglxs=1; nglys=1;        % 1x1 Gauss-Legendre quadrature for shear 
ngls=nglxs*nglys;        % number of sampling points per element for shear

%--------------------------------- % input data for nodal coordinate values %---------------------------------

gcoord=[0.0 0.0 0.0; 0.0 6.25 0.0; 0.0 12.5 0.0; 0.0 18.75 0.0;0.0 25.0 0.0; 
4.34 0.0 -0.38;4.34 6.25 -0.38;4.34 12.5 -0.38;4.34 18.75 -0.38;4.34 25.0 -0.38;     
8.55 0.0 -1.51;8.55 6.25 -1.51;8.55 12.5 -1.51;8.55 18.75 -1.51;8.55 25.0 -1.51;     
12.5 0.0 -3.35;12.5 6.25 -3.35;12.5 12.5 -3.35;12.5 18.75 -3.35;12.5 25.0 -3.35;     
16.1 0.0 -5.85;16.1 6.25 -5.85;16.1 12.5 -5.85;16.1 18.75 -5.85;16.1 25.0 -5.85];    
gcoord=12.0*gcoord;

%----------------------- % input data for nodal connectivity for each element %-----------------------

 nodes=[6 7 2 1; 7 8 3 2; 8 9 4 3; 9 10 5 4;
      11 12 7 6; 12 13 8 7; 13 14 9 8; 14 15 10 9;
      16 17 12 11; 17 18 13 12; 18 19 14 13; 19 20 15 14;
      21 22 17 16; 22 23 18 17; 23 24 19 18; 24 25 20 19];

%------------------------------------- % input data for boundary conditions %-------------------------------------

       bcdof=[1 3 5 6   7 11 12  13 17 18   19 23 24   25 26 28 29 30 ...
       31 33   56 58 60   61 63   86 88 90   91 93   116 118 120 ...
       121 122 123   146 148 150]; % constrained dofs
bcval=zeros(size(bcdof));     % whose described values are 0

%-------------------------------------- % initialization of matrices and vectors %------------------------------------

 ff=zeros(sdof,1);        % system force vector
 kk=zeros(sdof,sdof);     % system matrix
 disp=zeros(sdof,1);      % system displacement vector
 index=zeros(edof,1);     % index vector
 kinmtsb=zeros(3,edof);   % kinematic matrix for bending
 matmtsb=zeros(3,3);      % constitutive matrix for bending
 kinmtsm=zeros(3,edof);   % kinematic matrix for membrane
 matmtsm=zeros(3,3);      % constitutive matrix for membrane
 kinmtss=zeros(2,edof);   % kinematic matrix for shear
 matmtss=zeros(2,2);      % constitutive matrix for shear
 tr3d=zeros(edof,edof);   % transformation matrix

%---------------------------- % force vector %----------------------------

 ff(3)=-613.6;              % transverse force at node 1
 ff(9)=-1227.2;              % transverse force at node 2
 ff(15)=-1227.2;             % transverse force at node 3
 ff(21)=-1227.2;              % transverse force at node 4
 ff(27)=-613.6;              % transverse force at node 5
 ff(33)=-1227.2;             % transverse force at node 6
 ff(39)=-2454.4;              % transverse force at node 7
 ff(45)=-2454.4;              % transverse force at node 8
 ff(51)=-2454.4;             % transverse force at node 9
 ff(57)=-1227.2;              % transverse force at node 10
 ff(63)=-1227.2;              % transverse force at node 11
 ff(69)=-2454.4;             % transverse force at node 12
 ff(75)=-2454.4;              % transverse force at node 13
 ff(81)=-2454.4;              % transverse force at node 14
 ff(87)=-1227.2;             % transverse force at node 15
 ff(93)=-1227.2;              % transverse force at node 16
 ff(99)=-2454.4;              % transverse force at node 17
 ff(105)=-2454.4;             % transverse force at node 18
 ff(111)=-2454.4;              % transverse force at node 19
 ff(117)=-1227.2;             % transverse force at node 20
 (123)=-613.6;              % transverse force at node 21
 ff(129)=-1227.2;              % transverse force at node 22
 ff(135)=-1227.2;             % transverse force at node 23
 ff(141)=-1227.2;              % transverse force at node 24
 ff(147)=-613.6;             % transverse force at node 25

%-------------------- % computation of element matrices and vectors and their assembly %---------------------

% for bending and membrane stiffness

 [pointb,weightb]=feglqd2(nglxb,nglyb);     % sampling points & weights
 matmtsm=fematiso(1,emodule,poisson)*t;       % membrane material property
 matmtsb=fematiso(1,emodule,poisson)*t^3/12;  % bending material property
 [points,weights]=feglqd2(nglxs,nglys);     % sampling points & weights
 shearm=0.5*emodule/(1.0+poisson);          % shear modulus
 shcof=5/6;                                 % shear correction factor
 matmtss=shearm*shcof*t*[1 0; 0 1];         % shear material property

 for iel=1:nel           % loop for the total number of elements

 for i=1:nnel
 nd(i)=nodes(iel,i);         % extract connected node for (iel)-th element
 xcoord(i)=gcoord(nd(i),1);  % extract x value of the node
 ycoord(i)=gcoord(nd(i),2);  % extract y value of the node
 zcoord(i)=gcoord(nd(i),3);  % extract z value of the node
 end

 % 
 %  compute the local direction cosines and local axes
 %
 [tr3d,xprime,yprime]=fetransh(xcoord,ycoord,zcoord,nnel);
 %

 k=zeros(edof,edof);         % element matrix in local axes
 ke=zeros(edof,edof);        % element matrix in global axes
 km=zeros(edof,edof);        % element membrane matrix
 kb=zeros(edof,edof);        % element bending matrix
 ks=zeros(edof,edof);        % element shear matrix

 %------------------------------------------------------
 %  numerical integration for bending term
 %------------------------------------------------------

 for intx=1:nglxb
 x=pointb(intx,1);                  % sampling point in x-axis
 wtx=weightb(intx,1);               % weight in x-axis
 for inty=1:nglyb
 y=pointb(inty,2);                  % sampling point in y-axis
 wty=weightb(inty,2) ;              % weight in y-axis

 [shape,dhdr,dhds]=feisoq4(x,y);     % compute shape functions and
                                    % derivatives at sampling point

 jacob2=fejacob2(nnel,dhdr,dhds,xprime,yprime);  % compute Jacobian

 detjacob=det(jacob2);                 % determinant of Jacobian
 invjacob=inv(jacob2);                 % inverse of Jacobian matrix

 [dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob); % derivatives w.r.t.
                                               % physical coordinate

 kinmtsb=fekinesb(nnel,dhdx,dhdy);          % bending kinematic matrix
 kinmtsm=fekinesm(nnel,dhdx,dhdy);          % membrane kinematic matrix

 %--------------------------------------------
 %  compute bending element matrix
 %--------------------------------------------

 kb=kb+kinmtsb'*matmtsb*kinmtsb*wtx*wty*detjacob;
 km=km+kinmtsm'*matmtsm*kinmtsm*wtx*wty*detjacob;
 end
 end                   % end of numerical integration loop for bending term

 %------------------------------------------------------
 %  numerical integration for bending term
 %------------------------------------------------------

 for intx=1:nglxs
 x=points(intx,1);                  % sampling point in x-axis
 wtx=weights(intx,1);               % weight in x-axis
 for inty=1:nglys
 y=points(inty,2);                  % sampling point in y-axis
 wty=weights(inty,2) ;              % weight in y-axis

 [shape,dhdr,dhds]=feisoq4(x,y);     % compute shape functions and
                                    % derivatives at sampling point

 jacob2=fejacob2(nnel,dhdr,dhds,xprime,yprime);  % compute Jacobian

 detjacob=det(jacob2);                 % determinant of Jacobian
 invjacob=inv(jacob2);                 % inverse of Jacobian matrix

 [dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob); % derivatives w.r.t.
                                               % physical coordinate

 kinmtss=fekiness(nnel,dhdx,dhdy,shape);        % shear kinematic matrix

 %----------------------------------------
 %  compute shear element matrix
 %----------------------------------------

 ks=ks+kinmtss'*matmtss*kinmtss*wtx*wty*detjacob;   

 end
 end                      % end of numerical integration loop for shear term

 %--------------------------------
 %  compute element matrix
 %--------------------------------

 k=km+kb+ks;

 %-----------------------------------------------
 %  transform from local to global systems
 %-----------------------------------------------

 ke=tr3d'*k*tr3d;

 index=feeldof(nd,nnel,ndof);% extract system dofs associated with element

 kk=feasmbl1(kk,ke,index);  % assemble element matrices 

 end

 %-------------------------------------
 %   check the singular drilling dof
 %-------------------------------------

 for i=1:sdof
 if(abs(kk(i,i)) < 1e-5) 

 sum=0.0;
 for j=1:sdof
 sum=sum+abs(kk(i,j));
 end

 if (sum < 1e-5) 
 kk(i,i)=1;
 end

 end

 end

 %-----------------------------
 %   apply boundary conditions
 %-----------------------------

 [kk,ff]=feaplyc2(kk,ff,bcdof,bcval);

 %----------------------------
 %  solve the matrix equation
 %----------------------------

 disp=kk\ff;   

 num=1:1:sdof;
 displace=[num' disp]                  % print nodal displacements

 %--------------------------------------------------------------
 %                      Sub_routines
                     %%%%%%%%%%%%%%%%%

 function [point2,weight2]=feglqd2(nglx,ngly)

 %-------------------------------------------------------------------
 %  Purpose:
 %     determine the integration points and weighting coefficients
 %     of Gauss-Legendre quadrature for two-dimensional integration
 %
 %  Synopsis:
 %     [point2,weight2]=feglqd2(nglx,ngly) 
 %
 %  Variable Description:
 %     nglx - number of integration points in the x-axis
 %     ngly - number of integration points in the y-axis
 %     point2 - vector containing integration points   
 %     weight2 - vector containing weighting coefficients 
 %------------------------------------------------------------------- 

 %  determine the largest one between nglx and ngly

   if nglx > ngly
      ngl=nglx;
   else
      ngl=ngly;
   end

 %  initialization

   point2=zeros(ngl,2);
   weight2=zeros(ngl,2);

 %  find corresponding integration points and weights

[pointx,weightx]=feglqd1(nglx);     % quadrature rule for x-axis
[pointy,weighty]=feglqd1(ngly);     % quadrature rule for y-axis

%  quadrature for two-dimension

for intx=1:nglx                     % quadrature in x-axis
 point2(intx,1)=pointx(intx);
 weight2(intx,1)=weightx(intx);
end

 for inty=1:ngly                     % quadrature in y-axis
   point2(inty,2)=pointy(inty);
   weight2(inty,2)=weighty(inty);
  end
 _______________________________________________________________________
 function [matmtrx]=fematiso(iopt,elastic,poisson)

 %------------------------------------------------------------------------
 %  Purpose:
 %     determine the constitutive equation for isotropic material
 %
 %  Synopsis:
 %     [matmtrx]=fematiso(iopt,elastic,poisson) 
 %
 %  Variable Description:
 %     elastic - elastic modulus
 %     poisson - Poisson's ratio   
 %     iopt=1 - plane stress analysis
 %     iopt=2 - plane strain analysis
 %     iopt=3 - axisymmetric analysis
 %     iopt=4 - three dimensional analysis
 %------------------------------------------------------------------------

 if iopt==1        % plane stress
   matmtrx= elastic/(1-poisson*poisson)* ...
   [1  poisson 0; ...
   poisson  1  0; ...
   0  0  (1-poisson)/2];

elseif iopt==2     % plane strain
 matmtrx= elastic/((1+poisson)*(1-2*poisson))* ...
 [(1-poisson)  poisson 0; 
 poisson  (1-poisson)  0;
 0  0  (1-2*poisson)/2];

 elseif iopt==3     % axisymmetry
   matmtrx= elastic/((1+poisson)*(1-2*poisson))* ...
   [(1-poisson)  poisson  poisson  0; 
   poisson  (1-poisson)   poisson  0;
   poisson  poisson  (1-poisson)   0;
   0    0    0   (1-2*poisson)/2];

 else     % three-dimension
   matmtrx= elastic/((1+poisson)*(1-2*poisson))* ...
   [(1-poisson)  poisson  poisson   0   0    0; 
   poisson  (1-poisson)   poisson   0   0    0;
   poisson  poisson  (1-poisson)    0   0    0;
   0    0    0    (1-2*poisson)/2   0    0;
   0    0    0    0    (1-2*poisson)/2   0;
   0    0    0    0    0   (1-2*poisson)/2];

 end
______________________________________________________________________

 function [tr3d,xprime,yprime]=fetransh(xcoord,ycoord,zcoord,n)

 %--------------------------------------------------------------
 %  Purpose:
 %     Compute direction cosines between three-dimensional
 %    local and global coordinate axes
 %
 %  Synopsis:
 %     [tr3d,xprime,yprime]=fetransh(xcoord,ycoord,zcoord,n) 
 %
 %  Variable Description:
 %     xcoord - nodal x coordinates (4x1)   
 %     ycoord - nodal y coordinates (4X1)
 %     zcoord - nodal z coordinates (4X1)
 %     n - number of nodes per element
 %      tr3d - 3d transformation matrix from local to global axes
 %     xprime - coordinate in terms of the local axes (4x1)
 %     yprime - coordinate in terms of the global axes (4X1)
 %
 %  Note:
 %     The local x-axis is defined in the direction from the first node
 %     to the second node. Nodes 1, 2 and 4 define the local xy-plane.
 %     The local z-axis is defined normal to the local xy-plane.
 %     The local y-axis is defined normal to the x and z axes.
 %--------------------------------------------------------------------------

 %
 %  compute direction cosines
 %
 v12x=xcoord(2)-xcoord(1); 
 v12y=ycoord(2)-ycoord(1);
 v12z=zcoord(2)-zcoord(1);
 l12=sqrt(v12x^2+v12y^2+v12z^2);
 v23x=xcoord(3)-xcoord(2);
 v23y=ycoord(3)-ycoord(2);
 v23z=zcoord(3)-zcoord(2);
 l23=sqrt(v23x^2+v23y^2+v23z^2);
 v34x=xcoord(4)-xcoord(3);
 v34y=ycoord(4)-ycoord(3);
 v34z=zcoord(4)-zcoord(3);
 l34=sqrt(v34x^2+v34y^2+v34z^2);
 v14x=xcoord(4)-xcoord(1);
 v14y=ycoord(4)-ycoord(1);
 v14z=zcoord(4)-zcoord(1);
 l14=sqrt(v14x^2+v14y^2+v14z^2);
 v13x=xcoord(3)-xcoord(1); 
 v13y=ycoord(3)-ycoord(1);
 v13z=zcoord(3)-zcoord(1);
 l13=sqrt(v13x^2+v13y^2+v13z^2);
 v1tx=v12y*v14z-v12z*v14y;
 v1ty=v12z*v14x-v12x*v14z;
 v1tz=v12x*v14y-v12y*v14x;
 v1yx=v1ty*v12z-v1tz*v12y;
 v1yy=v1tz*v12x-v1tx*v12z;
 v1yz=v1tx*v12y-v1ty*v12x;
 vxx=v12x/l12;
 vxy=v12y/l12;
 vxz=v12z/l12;
 vyx=v1yx/sqrt(v1yx^2+v1yy^2+v1yz^2);
 vyy=v1yy/sqrt(v1yx^2+v1yy^2+v1yz^2);
 vyz=v1yz/sqrt(v1yx^2+v1yy^2+v1yz^2);
 vzx=v1tx/sqrt(v1tx^2+v1ty^2+v1tz^2);
 vzy=v1ty/sqrt(v1tx^2+v1ty^2+v1tz^2);
 vzz=v1tz/sqrt(v1tx^2+v1ty^2+v1tz^2);

 %
 %   transformation matrix
 % 
 for i=1:2*n
 i1=(i-1)*3+1;
 i2=i1+1;
 i3=i2+1;
 tr3d(i1,i1)=vxx;
 tr3d(i1,i2)=vxy;
 tr3d(i1,i3)=vxz;
 tr3d(i2,i1)=vyx;
 tr3d(i2,i2)=vyy;
 tr3d(i2,i3)=vyz;
 tr3d(i3,i1)=vzx;
 tr3d(i3,i2)=vzy;
 tr3d(i3,i3)=vzz;
 end

 %   
 %  compute nodal values in terms of local axes
 %
 alpa213=acos((l12^2+l13^2-l23^2)/(2*l12*l13));
 alpa314=acos((l13^2+l14^2-l34^2)/(2*l13*l14));
 alpa41y=2*atan(1)-alpa213-alpa314;
 xprime(1)=0;  yprime(1)=0;
 xprime(2)=l12;  yprime(2)=0;
 xprime(3)=l13*cos(alpa213);  yprime(3)=l13*sin(alpa213);
 xprime(4)=l14*sin(alpa41y);  yprime(4)=l14*cos(alpa41y);
_____________________________________________________________________

 function [shapeq4,dhdrq4,dhdsq4]=feisoq4(rvalue,svalue)

 %------------------------------------------------------------------------
 %  Purpose:
 %     compute isoparametric four-node quadilateral shape functions
 %     and their derivatves at the selected (integration) point
 %     in terms of the natural coordinate 
 %
 %  Synopsis:
 %     [shapeq4,dhdrq4,dhdsq4]=feisoq4(rvalue,svalue)  
 %
 %  Variable Description:
 %     shapeq4 - shape functions for four-node element
 %     dhdrq4 - derivatives of the shape functions w.r.t. r
%     dhdsq4 - derivatives of the shape functions w.r.t. s
%     rvalue - r coordinate value of the selected point   
%     svalue - s coordinate value of the selected point
%
%  Notes:
%     1st node at (-1,-1), 2nd node at (1,-1)
%     3rd node at (1,1), 4th node at (-1,1)
%------------------------------------------------------------------------

 % shape functions

 shapeq4(1)=0.25*(1-rvalue)*(1-svalue);
 shapeq4(2)=0.25*(1+rvalue)*(1-svalue);
 shapeq4(3)=0.25*(1+rvalue)*(1+svalue);
 shapeq4(4)=0.25*(1-rvalue)*(1+svalue);

 % derivatives

 dhdrq4(1)=-0.25*(1-svalue);
 dhdrq4(2)=0.25*(1-svalue);
 dhdrq4(3)=0.25*(1+svalue);
 dhdrq4(4)=-0.25*(1+svalue);

 dhdsq4(1)=-0.25*(1-rvalue);
 dhdsq4(2)=-0.25*(1+rvalue);
 dhdsq4(3)=0.25*(1+rvalue);
 dhdsq4(4)=0.25*(1-rvalue);
________________________________________________________________________

 function [jacob2]=fejacob2(nnel,dhdr,dhds,xcoord,ycoord)

 %------------------------------------------------------------------------
 %  Purpose:
 %     determine the Jacobian for two-dimensional mapping
 %
 %  Synopsis:
 %     [jacob2]=fejacob2(nnel,dhdr,dhds,xcoord,ycoord) 
 %
 %  Variable Description:
 %     jacob2 - Jacobian for one-dimension
 %     nnel - number of nodes per element   
 %     dhdr - derivative of shape functions w.r.t. natural coordinate r
 %     dhds - derivative of shape functions w.r.t. natural coordinate s
 %     xcoord - x axis coordinate values of nodes
 %     ycoord - y axis coordinate values of nodes
 %------------------------------------------------------------------------

jacob2=zeros(2,2);

 for i=1:nnel
 jacob2(1,1)=jacob2(1,1)+dhdr(i)*xcoord(i);
 jacob2(1,2)=jacob2(1,2)+dhdr(i)*ycoord(i);
 jacob2(2,1)=jacob2(2,1)+dhds(i)*xcoord(i);
 jacob2(2,2)=jacob2(2,2)+dhds(i)*ycoord(i);
 end
________________________________________________________________________

 function [dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob)

 %------------------------------------------------------------------------
 %  Purpose:
 %     determine derivatives of 2-D isoparametric shape functions with 
 %     respect to physical coordinate system
 %
 %  Synopsis:
 %     [dhdx,dhdy]=federiv2(nnel,dhdr,dhds,invjacob)  
 %
 %  Variable Description:
 %     dhdx - derivative of shape function w.r.t. physical coordinate x
 %     dhdy - derivative of shape function w.r.t. physical coordinate y
 %     nnel - number of nodes per element   
 %     dhdr - derivative of shape functions w.r.t. natural coordinate r
 %     dhds - derivative of shape functions w.r.t. natural coordinate s
 %     invjacob - inverse of 2-D Jacobian matrix
 %------------------------------------------------------------------------

 for i=1:nnel
 dhdx(i)=invjacob(1,1)*dhdr(i)+invjacob(1,2)*dhds(i);
 dhdy(i)=invjacob(2,1)*dhdr(i)+invjacob(2,2)*dhds(i);
 end
_________________________________________________________________________
function [kinmtsb]=fekinesb(nnel,dhdx,dhdy)

 %--------------------------------------------------------------------------
  %  Purpose:
 %     determine the kinematic matrix expression relating bending curvatures 
 %     to rotations and displacements for shear deformable plate bending
 %
 %  Synopsis:
 %     [kinmtsb]=fekinesb(nnel,dhdx,dhdy) 
 %
 %  Variable Description:
 %     nnel - number of nodes per element
 %     dhdx - derivatives of shape functions with respect to x   
 %     dhdy - derivatives of shape functions with respect to y
 %--------------------------------------------------------------------------

 for i=1:nnel
 i1=(i-1)*6+1;  
 i2=i1+1;
 i3=i2+1;
 i4=i3+1;
 i5=i4+1;
 i6=i5+1;
 kinmtsb(1,i5)=dhdx(i);
 kinmtsb(2,i4)=-dhdy(i);
 kinmtsb(3,i5)=dhdy(i);
 kinmtsb(3,i4)=-dhdx(i);
 kinmtsb(3,i6)=0;
 end
_________________________________________________________________________
function [kinmtsm]=fekinesm(nnel,dhdx,dhdy)

 %------------------------------------------------------------------------
 %  Purpose:
  %     determine the kinematic equation between strains and displacements
 %     for two-dimensional solids
 %
 %  Synopsis:
 %     [kinmtsm]=fekinesm(nnel,dhdx,dhdy) 
 %
 %  Variable Description:
 %     nnel - number of nodes per element
 %     dhdx - derivatives of shape functions with respect to x   
 %     dhdy - derivatives of shape functions with respect to y
 %------------------------------------------------------------------------

 for i=1:nnel
 i1=(i-1)*6+1;  
 i2=i1+1;
 i3=i2+1;
 i4=i3+1;
 i5=i4+1;
 i6=i5+1;
 kinmtsm(1,i1)=dhdx(i);
 kinmtsm(2,i2)=dhdy(i);
 kinmtsm(3,i1)=dhdy(i);
 kinmtsm(3,i2)=dhdx(i);
 kinmtsm(3,i6)=0.0;
 end
__________________________________________________________________________
 function [kinmtss]=fekiness(nnel,dhdx,dhdy,shape)

 %------------------------------------------------------------------------
 %  Purpose:
 %     determine the kinematic matrix expression relating shear strains 
 %     to rotations and displacements for shear deformable plate bending
 %
 %  Synopsis:
 %     [kinmtss]=fekiness(nnel,dhdx,dhdy,shape) 
 %
 %  Variable Description:
 %     nnel - number of nodes per element
 %     dhdx - derivatives of shape functions with respect to x   
 %     dhdy - derivatives of shape functions with respect to y
 %     shape - shape function
 %------------------------------------------------------------------------

 for i=1:nnel
 i1=(i-1)*6+1;  
 i2=i1+1;
 i3=i2+1;
 i4=i3+1;
 i5=i4+1;
 i6=i5+1;
 kinmtss(1,i3)=dhdx(i);
 kinmtss(1,i5)=shape(i);
 kinmtss(2,i3)=dhdy(i);
 kinmtss(2,i4)=-shape(i);
 kinmtss(2,i6)=0.0;
 end
_________________________________________________________________________
function [index]=feeldof(nd,nnel,ndof)
%----------------------------------------------------------
%  Purpose:
%     Compute system dofs associated with each element 
%
%  Synopsis:
%     [index]=feeldof(nd,nnel,ndof)
%
%  Variable Description:
%     index - system dof vector associated with element "iel"
%     iel - element number whose system dofs are to be determined
%     nnel - number of nodes per element
%     ndof - number of dofs per node 
%-----------------------------------------------------------

 edof = nnel*ndof;
   k=0;
   for i=1:nnel
     start = (nd(i)-1)*ndof;
       for j=1:ndof
         k=k+1;
         index(k)=start+j;
       end
   end

 ________________________________________________________________________
   function [kk]=feasmbl1(kk,k,index)
 %----------------------------------------------------------
 %  Purpose:
 %     Assembly of element matrices into the system matrix
 %
 %  Synopsis:
 %     [kk]=feasmbl1(kk,k,index)
 %
 %  Variable Description:
 %     kk - system matrix
 %     k  - element matri
 %     index - d.o.f. vector associated with an element
 %-----------------------------------------------------------

 edof = length(index);
 for i=1:edof
   ii=index(i);
     for j=1:edof
       jj=index(j);
         kk(ii,jj)=kk(ii,jj)+k(i,j);
     end
 end

________________________________________________________________________

 function [kk,ff]=feaplyc2(kk,ff,bcdof,bcval)

 %----------------------------------------------------------
 %  Purpose:
 %     Apply constraints to matrix equation [kk]{x}={ff}
 %
 %  Synopsis:
 %     [kk,ff]=feaplybc(kk,ff,bcdof,bcval)
 %
 %  Variable Description:
 %     kk - system matrix before applying constraints 
 %     ff - system vector before applying constraints
 %     bcdof - a vector containging constrained d.o.f
 %     bcval - a vector containing contained value 
 %
 %     For example, there are constraints at d.o.f=2 and 10
 %     and their constrained values are 0.0 and 2.5, 
 %     respectively.  Then, bcdof(1)=2 and bcdof(2)=10; and
 %     bcval(1)=1.0 and bcval(2)=2.5.
 %-----------------------------------------------------------

 n=length(bcdof);
 sdof=size(kk);

 for i=1:n
    c=bcdof(i);
    for j=1:sdof
       kk(c,j)=0;
    end

    kk(c,c)=1;
    ff(c)=bcval(i);
 end
________________________________________________________________________

[Walter Roberson]

Please read the guide to tags and retag this question. See http://www.mathworks.co.uk/matlabcentral/answers/43073-a-guide-to-tags

Answer 1

[Walter Roberson]

The error was made on line 71 by Colonel Mustard with the Wrench.

No? How about a Clue ?

[Deviprakash Upadhyaya]

Sorry for the problem created.I have edited the question and made it proper.Will be getting a warning message in " disp=kk\ff; ".Problem with inversing the system matrix " kk ".