$0.1$ is not exactly representable in double precision arithmetic, and it will be rounded. The double precision representation of $0.1$ is
$\verb|00111111 10111001 10011001 10011001 10011001 10011001 10011001 10011010|,$
which is $0.100000000000000005551115123126$. As you can see, if you compute $h$ by subtracting "$0.1$" from it repeatedly, you ...
To make my comments more explicit:
Try with the following code for the DDE function that eliminates P
function dydt = dde(t,y,yd,A)
n = length(A);
dydt = diag(A).*y;
k = 0;
for i = 1:n
for j = (i+1):n
dydt(i) += A(i,j)*yd(j,k);
dydt(j) += A(i,j)*yd(i,k);
k = k+1;
This eliminates the $O(n^3)$ matrix-vector ...
Let's take it from the beginning. I was unable to run your code as it must be missing some input, but I see that your input to the differential equations is the current $I$ which you control as a user, and it looks something like this:
You will notice that this input is discontinuous, and thus non-differentiable. Unfortunately matlab does not have any ...
The elemental stiffness matrix must be always singular because while deriving it we do not impose any constraints or boundary conditions. Thus inverting the stiffness matrix to solve for displacements/position-vectors/degrees-of-freedom should yield indeterminate results.
I think that what you are looking for is a Computer Algebra System (CAS). You could find a list in Wikipedia.
I don't know which one suits your particular interest, although Singular seems to be a good candidate.