Hi Sahil,
This can be done by for loop but it is maybe better to do the whole thing with index manipulation. In the written-out expressions for dy, looking down the columns of indices there are some circular shifts between some of the columns. There are two kinds of indices, those that increment by 1, and the y indices that increment by 3. Denoting the first set of indices with n, then
n1 = (1:20)';
n2 = circshift(n1,-1);
n61 = n1 +60;
n62 = circshift(n61,-1);
n80 = circshift(n61,1);
% as a check, compare the indices for o( ) and sin(y( )) with the written-out version
nset = [n1 n2 n61 n62 n80]
nset =
1 2 61 62 80
2 3 62 63 61
3 4 63 64 62
4 5 64 65 63
5 6 65 66 64
6 7 66 67 65
7 8 67 68 66
8 9 68 69 67
9 10 69 70 68
10 11 70 71 69
11 12 71 72 70
12 13 72 73 71
13 14 73 74 72
14 15 74 75 73
15 16 75 76 74
16 17 76 77 75
17 18 77 78 76
18 19 78 79 77
19 20 79 80 78
20 1 80 61 79
Here n80 is called that because its value is 80 in the first row, and similarly for the others, first row.
For the increment-by-3 indices,
j2 = 3*(1:20)-1;
j5 = circshift(j2,-1);
j8 = circshift(j2,-2);
j59 = circshift(j2,1);
% as a check, compare these y( ) indices with the written-out version
jset = [j2; j5; j8; j59]'
jset =
2 5 8 59
5 8 11 2
8 11 14 5
11 14 17 8
14 17 20 11
17 20 23 14
20 23 26 17
23 26 29 20
26 29 32 23
29 32 35 26
32 35 38 29
35 38 41 32
38 41 44 35
41 44 47 38
44 47 50 41
47 50 53 44
50 53 56 47
53 56 59 50
56 59 2 53
59 2 5 56
NOTE that your expression for dy(80) contains two values of 57 which by all rights look like they should be 59.
Given the indices, all of which are vectors of length 20, you can calculate all of the dy in one go:
dy(n61) = o(1,n2) - o(1,n1) ...
- (k./tp).*(y(j2) ./y(j5)).*sin(y(n61)) - (k./tp).*(y( j5) ./y(j2)).*sin(y(n61)) ...
+ (k./tp).*(y(j8) ./y(j5)).*sin(y(n62)) + (k./tp).*(y(j59) ./y(j2)).*sin(y(n80));
for the plotting:
figure(1)
plot(T./tp,(Y(:,61)));
hold on
for m = 62:80
plot(T./tp,(Y(:,m)));
end
hold off
