Problème avec matlab

Bonjour à tous
dans le cadre de mon travail de diplôme je me suis heurté à un problème de matériel et j'ai dû contourner l'obstacle en modélisant une plaque et e seul langage abordable pour moi est Matlab que je ne maîtrise pas bien mais j'évolue à un rythme pas très rapide c'est la raison pour laquelle je suis là
J'ai découvert un projet sur le net qui répond à ma demande j'ai donc commencé à résoudre ce projet que je vous mets en attaché.et ce que j'ai fait je ne comprends pas pourquoi je n'ai pas tous les modes affichés et je bloque pour la suite avec le dirac.
Merci pour votre aide elle m'est précieuse.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
h=0.02; %épaisseur de la plaque
a=0.8; %Longueur de la plaque
b=0.8; %Largeur de la plaque
rho=540; %Masse volumique de la plaque
E=11*10^9; %Module de Young
khi=0.5; %Coeficient de Poisson
sigmap=0.003; %l'amortissement
m=rho*h; %masse surfacique
D=(E*h^3)/(12*(1-khi^2)); %la raideur en flexion de la plaque
p=0:5;
q=0:5;
omegapq=pi^2*sqrt(D/m)*((p.*p/a^2)+(q.*q/b^2)); %Fréquence des oscillations du mode considéré
sort(omegapq);
p=2;
q=1;
x=0:0.01:1;
y=0:0.01:1;
phixy=2/sqrt(a*b)*sin(p*pi*x/a)'*sin(pi*q*y/b); % Fonctions propres
mesh(x,y,phixy);
meshc(x,y,phixy);
repmat(phixy,10,10);
%omegapp=1.0e+004*0.0928;
kf=sqrt(sqrt((rho*h*(omegapq.*omegapq))/D)) % Le nombre d'onde de flexion de la plaque couplée au fluide
plot(kf)
plot(omegapq)
% Dans ce qui suit nous essayons de reconstruire le signal en affichant les
% différents modes et en les sommant pour avoir la figure finale du
% déplacement de l'onde
phixy = []
figure
i = 1
for p = 0:5
    for q = 0:5
        if p == 0 && q == 0, 
            phixy = 2/sqrt(a*b)*sin(p*pi*x/a)'*sin(pi*q*y/b);
%             subplot(5,5,i),mesh(x,y,phixy)
%             i=i+1;
        else
            subplot(5,5,i),mesh(x,y,2/sqrt(a*b)*sin(p*pi*x/a)'*sin(pi*q*y/b))
            title(['(p,q) = (',num2str(p),',',num2str(q),')'])
            phixy = phixy + 2/sqrt(a*b)*sin(p*pi*x/a)'*sin(pi*q*y/b);
            if i <25,i=i+1;end
        end
    end
end
figure,mesh(x,y,phixy) % La somme affichée.
omegapv=omegapq*sqrt(1-sigmap^2);
q = 2;
p=1:5;
t=0:5;
x=.34;
y=.20;
Fp = 1/m* 2/sqrt(a*b)*sin(p*pi*x/a)'*sin(pi*q*y/b)*dirac(t);