Functions that don't appear in the ODE system can also be plotted with odeplot. So, you can do

plots:-odeplot(Sol, [P,V](x), x= 0..10)

even though only P(x) appears in the differential equation(s). So, there's no need to change your original ODE(s).

The above command plots one curve, V vs P, not two separate curves.

Just to summarize the essential point, which may be difficult to see amidst the dense numeric code: The parameter of animate is inherently real-valued; it must be somehow converted to an integer before being used as an index to your table SOLNSuy. The alternative to animate is to create an indexable container of static plots and pass it to display(..., insequence).

What you found, unapply, is the most widely recommended way to do it.

The solution is simply the obvious: f(x/a).

I don't understand exactly why what you had doesn't work, but I strongly prefer using eval rather than assign. By simply making that change, it now works. Replace your last line of code with

Sol:= {};
for i from 0 to N do   
    Sol:= Sol union {dsolve({cond[1][i], eval(equ[2][i], Sol)}, {f[i](x)})}
end do;

It can be done like this:

C:= lcoeff(y^2 - x^2/y, order= plex(x,y), 't'):
C*t;
                       -x^2/y

In many top-level Maple commands, the concept of polynomial is generalized a bit to allow for negative exponents. 

You're making two mistakes: The first is mixing up integer and polynomial computation; the second is mixing up exact and floating-point computation. 

Integer vs polynomial: You're making exactly the same mistake as you were making regarding gcd. The mod commands (there's a whole family of them) are intended primarily for polynomials and matrices over finite rings (including finite fields). The command for computing the remainder of integer division is irem. Do you see a pattern here? Top-level commands for pure-integer arithmetic begin with i. Your command x mod 3 returns x because x is a polynomial; it's not a placeholder for an integer.

Exact vs floating-point: Neither of these exact commands is appropriate for use as the first argument to plot, which is necessarily based on floating-point computation. 

Replace sum with add. The problem with sum is that the derivative of u[i] is evaluated before i has a value. Since u[i] hasn't been defined in terms of x for generic i, its derivative is 0. On the other hand, add evaluates its index variable to a value in the specified range before evaluating its expression. 

One is tricked into believing that the first assembly in VV's Answer is a right triangle, making the area 32-1/2. But actually that assembly is a quadrilateral. The slope of the hypotenuse of the green triangle is 2/5, while that of the red triangle is 3/8. The difference of the angles is only about 1.2 degrees, or a relative difference of about 6%. Thus the supposed hypotenuse of the overall triangle is actually 2 sides of a quadrilateral, and its area is actually 32.

This is not a bug. The 2 represents the precision at which the computations are done, not the precision of the final result. The number 1014.1 can't be represented adequately in 2 digits.

Here's how to do it.
 

Download 3dModelFit.mw

The following is a variation of VV's procedure that uses the low-level operations of LinearAlgebra:-Modular for efficiency. I believe that it'll run about 2-1/2 times faster---not a huge improvement, but perhaps worthwhile. Please test it on a matrix known to be MDS.

IsMDS:= proc(
    A::Matrix(square),
    r::And(posint, satisfies(r-> irem(upperbound(A,1), r) = 0))
)
uses LAM= LinearAlgebra:-Modular;
local
    n:= iquo(upperbound(A,1), r), 
    rng:= t-> (t-1)*r+1..t*r,
    M:= LAM:-Mod(2, A, float[8]),
    k, kr, P, ip, jp, det
;
    for k to n do
        kr:= k*r;
        P:= combinat:-choose(n,k);
        for ip in P do 
            for jp in P do
                LAM:-RowReduce(
                    2, M[rng~(ip), rng~(jp)], kr$3, 'det', 0$4, false
                );
                if det=0 then return false fi
            od
        od
    od;
    true
end proc
:    

My try:

restart:
d:= proc(m::integer)
option remember;
local i;
    -expand(add(xi[i]*thisproc(m-i+1), i= 2..m))
end proc:
d(1):= 1:

DD:= m-> rtable((1..m)$2, (i,j)-> d(j-i+1), subtype= Matrix): 
DD(4);

The command gcd is meant for polynomials. As polynomials, a and b have no common factors, so their GCD is 1. The command for the GCD of integers is igcd.

You make a, b, and x the parameters of a procedure and pass to them two expressions and a variable name. To do RK4, we also need an initial condition, a final value of x, and a number of steps. Like this:

MyRK4:= proc(
    a::algebraic, b::algebraic, x::name,
    x0::realcons, y0::realcons, xf::realcons, n::posint
)
local y, h:= (xf-x0)/n, k;
    dsolve(
        {diff(y(x),x) = 10*(y(x)-a)/b, y(x0) = y0}, numeric,
        'method'= 'classical'['rk4'], 'stepsize'= h,
        'output'= Array(1..n, k-> x0+k*h)
    )[2][1]
end proc
:
MyRK4(ln(x), x, x, 1, 1, 2, 10);