lor1.m

四面体有限元程序,单元刚度的计算及总刚度矩征集成

           function varargout = lor1(t,y,flag)
           if nargin == 0,    
              flag = 'demo';  
           end

           switch flag
           case ''                        % Return dy/dt = f(t,y).
              varargout{1} = f(t,y);

           case 'demo'                    % Run a demo.
              demo;
           otherwise
              error(['Unknown flag ''' flag '''.']);
           end

           function dydt = f(t,y)
           A = [ -8/3     0     y(2)  
                    0   -10      10 
                 -y(2)   35      -1  ];
           dydt = A*y;

           function demo
           axis([10 50 -50 50 -50 50]);
           view(3)
           hold on
           title('Lonrenz Attractor')

           options =odeset('OutputFcn','odephas3') ;              
           [t,y]=ode23('lor1',[0,20],[0,0,eps],options);