I am trying to use the MATLAB Damped Cantilever Beam example with my own parameters but it is not working. I am getting an error Error using pde.PDERes​ults>class​ifyPDEMode​l (line 195) Solution does not correspond to the provided model.

Here is how my modified code(taken from CantileverBeamDampedTransient.m example file) looks like:

    %%Beam Dimensions and Material Properties
    % The code uses Imperial units.
    % The beam is 188.97 inches long and 1.18 inches thick.
    width = 188.97; 
    height = 1.18;

The material is aluminium.

    E = 9.86e6; 
    nu = 0.36; 
    rho = 0.0975;

Calculate the coefficient matrix from the material properties.

    G = E/(2.*(1+nu));
    mu = 2.0*G*nu/(1-nu);

    %%Lowest Vibration Frequency From Beam Theory
    I = height^3/12; 
    A = height;

From beam theory, there is a simple expression for the lowest vibration frequency of the cantilever beam.

    eigValAnalytical = 1.8751^4*E*I/(rho*A*width^4); 
    freqAnalytical = sqrt(eigValAnalytical)/(2*pi);

    %%Create a PDE Analysis Model
    % Create a PDEModel with two independent variables to represent the
    % analysis.
    numberOfPDE = 2;
    model = createpde(numberOfPDE);

    %%Create the Geometry
    % Create a simple rectangular geometry.
    gdm = [3;4;0;width;width;0;0;0;height;height];
    g = decsg(gdm,'S1',('S1')');

Provide the model with a definition of the geometry.

    geometryFromEdges(model,g);

    %%Equation Coefficients
    % The equation coefficients are derived from the material properties.
    c = [2*G+mu;0;G;0;G;mu;0;G;0;2*G+mu];
    f = [0;0];
    a = 0;
    m = rho;
    specifyCoefficients(model,'m',m,'d',0,'c',c,'a',a,'f',f);
    %%Boundary Conditions
    % Specify the following boundary condition to clamp (displacements equal
    % zero) the left beam-edge.
    applyBoundaryCondition(model,'dirichlet','Edge',4,'u',[0 0]);

    %%Mesh Generation
    % Define a maximum element size (5 elements through the beam thickness).
    hmax = height/5;
    msh=generateMesh(model,'Hmax',hmax,'MesherVersion','R2013a');

    %%Calculation of Vibration Modes and Frequencies
    % Use |solvepdeeig| and then compute the lowest-frequency vibration mode.
    res = solvepdeeig(model, [0,1e6]');
    eigenVec = res.Eigenvectors;
    eigenVal = res.Eigenvalues;
    uEigenTip = eigenVec(2,2);
    u0TipDisp = .1;

    u0 = u0TipDisp/uEigenTip*eigenVec;

Color plot of y-displacement of the lowest-frequency vibration mode.

    freqNumerical = sqrt(eigenVal(1))./(2*pi);
    longestPeriod = 1/freqNumerical;
    % Calculate the solution for three full periods.
    tlist = 0:longestPeriod/100:5*longestPeriod;
    % Create a function handle that can be used to provide initial conditions.
    R = createPDEResults(model, eigenVec(:));
    ice = icEvaluator(R);

Now solve the system with damping equal to 3% of critical for this lowest vibration frequency.

    fem = assembleFEMatrices(model);
    zeta = .03; 
    omega = 2*pi*freqNumerical; 
    alpha = 2*zeta*omega;
    dampmat = alpha*fem.M;
    %%Transient Analysis, Initial Displacement From Static Solution
    %
    % It would be very unusual for a structure to be loaded such that the
    % displacement is equal to a multiple of one of its vibration modes.
    % In this more realistic example, we solve the transient equations
    % with the initial displacement shape calculated from the static solution
    % of the cantilever beam with a vertical load at the tip.
    %
    % Perform a static analysis with vertical tip load equal one.
    % Follow the model building steps as before.
    P = 134.885;
    pdeTipLoad = createpde(2);
    pg = geometryFromEdges(pdeTipLoad,g);

Specify the equation to solve.

    specifyCoefficients(pdeTipLoad,'m',0,'d',0,'c',c,'a',a,'f',f);
    tipLoad = applyBoundaryCondition(pdeTipLoad,'neumann','Edge',2,'g',[0 P/height]);
    clampedEdge = applyBoundaryCondition(pdeTipLoad,'dirichlet','Edge',4,'u',[0,0]);
    msh=generateMesh(pdeTipLoad,'Hmax',hmax,'MesherVersion','R2013a');

    statres = solvepde(pdeTipLoad);

To make comparison with the eigenvector case clearer, we will also scale this static solution so that the maximum y-displacement at the tip equals .1 inches.

    u = statres.NodalSolution;
    uEigenTip = u(2,2); 
    u0TipDisp = 0.1;
    u0 = u0TipDisp/uEigenTip*u;

Calculate the undamped solution with the initial displacement from the static analysis.

    specifyCoefficients(model, 'm', m, 'd', 0, 'c', c, 'a', a, 'f', f);

Set the initial conditions and solve the system.

    R = createPDEResults(model, u0(:));
    ice = icEvaluator(R);
    setInitialConditions(model, @ice.computeIC, 0);
    tres = solvepde(model,tlist);

The displacement is no longer a pure sinusoidal wave. The static solution that we are using as the initial conditions is similar to the lowest eigenvector but higher-frequency modes are also contributing to the transient solution.

    titl = 'Initial Displacements from Static Solution, Undamped Solution';
    cantileverBeamTransientPlot(tres,titl);

Calculate the damped solution with the initial displacement from the static analysis. The damping matrix is the same as that used in the eigenvector case, above.

    specifyCoefficients(model, 'm', m, 'd', dampmat, 'c', c, 'a', a, 'f', f);
    tres = solvepde(model,tlist);

Plot the tip displacement from this solution as a function of time. Again we superimpose a curve of the damped amplitude as a function of time obtained from an analytical solution to the single degree-of-freedom ODE. Because the initial condition differs from the lowest eigenvector, this analytical solution only approximates the amplitude of the PDE solution.

    titl = 'Initial Displacements from Static Solution, Damped Solution';
    cantileverBeamTransientPlot(tres,titl);
    %hold on;
    %plot(tlist,u.*exp(-zeta*omega*tlist),'Color','r');
    %legend('PDE','ODE Amplitude Decay','Location','southeast');
    %%Utility Plot Function
    % Utility function for creating the plots of tip y-displacement as a 
    % function of time.
    type cantileverBeamTransientPlot.m

    %%Function Handle for Specifying Initial Condition (IC)
    % The evaluation function returns the value of the IC at any point within
    % the mesh. A results object represents the values and the interpolation
    % function that it provided is used to compute the ICs.
    type icEvaluator.m