restart;
a:=(x,v)->-x-0.1*v; #Example only
xo := 1; vo := 0;
sys:={diff(x(t), t) = v(t), diff(v(t), t) = a(x(t), v(t)),x(0)=xo,v(0)=vo}; #Your system
dsolve(sys); #For comparison later
X:=subs(%,x(t));
plot(X,t=0..50);
#Now the corrected algorithm
xo := 1; vo := 0:
A := Matrix():
dt := 0.1e-2:
A(1, 1) := xo; A(1, 2) := vo; A(1, 3) := 0;
for i from 2 to 1001 do
   xn := xo+vo*dt; vn := vo+a(xo,vo)*dt;
   dx1 := vn*dt; dv1 := a(xo,vo)*dt;
   dx2 := (vn+(1/2)*dv1)*dt; dv2 := a(xn+(1/2)*dx1,vn+(1/2)*dv1)*dt;
   dx3 := (vn+(1/2)*dv2)*dt; dv3 := a(xn+(1/2)*dx2, vn+(1/2)*dv2)*dt;
   dx4 := (vn+dv3)*dt; dv4 := a(xn+dx3, vn+dv3)*dt;
   xn := evalf(xo+(dx1+2*dx2+2*dx3+dx4)*(1/6)); vn := evalf(vo+(dv1+2*dv2+2*dv3+dv4)*(1/6));
   A(i, 1) := xn; A(i, 2) := vn; A(i, 3) := ( i-1)*dt;
   xo := xn;
   vo := vn
end do:
A[1,..];
#The submatrix consisting of the third and first column can be obtained in two ways:
#LinearAlgebra:-SubMatrix(A,[1..-1],[3,1]);
# or
#A[1..-1,[3,1]];
plot(A[1..-1,[3,1]]);
prk:=%:
plot(X,t=0..1,color=blue):
plots:-display(prk,%);
#Comparing values at t = 1:
A[-1,..];
evalf(eval(X,t=1));

Edited: 1000 changed to 1001 and A(i, 3) :=  i*dt changed to A(i, 3) := ( i-1)*dt

Integrating by parts twice we obtain:

Int(cos(t^3)/(t+2), t = -infinity .. x) = (1/3)*sin(x^3)/((x+2)*x^2)+O(1/x^6)

as is seen here:

restart;
with(IntegrationTools):
A:=Int(cos(t^3)/(t+2), t = -infinity .. x);
#By changing variable by hand we get:
B:=Int(cos(t^3)/(2-t), t = -x..infinity);
C:=Change(B,t=s^(1/3));
subs(cos(s)=1,GetIntegrand(C));
Parts(C,%); #First integration by parts
res:=map(simplify,%) assuming x < -2;
C2:=op(2,res);
subs(sin(s)=1,GetIntegrand(C2));
Parts(C2,%); #Second integration by parts
res2:=map(simplify,%) assuming x < -2;
#The integral is bounded by a constant times:
Int(1/(s-8)^3, s = -x^3 .. infinity);
value(%) assuming x < -2;
#Thus we get
A=op(1,res)+O(x^(-6));

#And clearly you can continue like that. With just one more integration by parts we can see that
Int(cos(t^3)/(t+2), t = -infinity .. x) =
(1/3)*sin(x^3)/((x+2)*x^2)-(1/9)*(3*x+4)*cos(x^3)/((x+2)^2*x^5)+O(1/x^9)

Your system seems to be {diff(x(t), t) = v(t), and diff(v(t), t) = a(x(t), v(t))}, where a is a function of two variables. However, that function is never made concrete.

The multiplication sign in Maple is *, not a dot (used for matrix multiplication though).

What are v and a in the lines xn := xo+v.dt; and vn := vo+a.dt; ?

There may be more problems.

The multiplication sign in Maple is *. Replace all dots by * :

restart;
e1 := 0 = a*Pyrcyt-a*Pyrmt+(b-c)*Malmt;
e2 := 0 = d*Malcyt+e*Citmt-(d+b+c)*Malmt;
e3 := 0 = f*Citcyt+c*Malmt-(e+f)*Citmt;
e4 := 0 = 2*e*Citmt+(b+c)*Malmt-Resp;
eqns := {e1, e2, e3, e4};
solve(eqns, {Citcyt, Citmt, Malmt, Pyrmt});

There are infinitely many solutions.

For each constant a<>0 the function

f:=x->sin(Pi/2/a*x);
#is a solution:
expand(f(x+y));
simplify(f(x)*f(a-y)+f(y)*f(a-x));
f(2013^(1/2013));

(These solutions were just found from recalling the well known addition formula for the sine function).

And by the way: By observation the constant function f(x) = 1/2 is a solution too (besides the trivial solution f(x)=0).

Use add instead of sum and use a procedure in fsolve:
fsolve(n->add(length(i), i = 1 .. n)- 2893, 500..2000);

The output (in my case anyway) was 1000.862939, but notice that

add(g(i),i=1..3.4567) = add(g(i),i=1..3);

Something like this works:

L:=seq(1/n,n=1..100):
f:=(M,N)->M+N+M*N;
while nops([L])>1 do
   n,m:=op(combinat:-randcomb(nops([L]),2)); #Two random indices in the range 1..nops([L])
   L:=f(L[n],L[m]),L[1..n-1],L[n+1..m-1],L[m+1..-1]
end do:
L;

Automatic simplification turns  (5+3*x)/15 into  1/3 + x/5.
Try

'(5+3*x)/15';

Thus any conceivable simplification procedure used on 1/3 + x/5 would fail to produce what you want unless it produces output like e.g.

``(5+3*x)/15;
(5+3*x)/``(15);

You have two equations and three unknowns.

Thus given e.g. kappa you can solve for alpha and beta using fsolve:

fsolve(eval({f,g},kappa=1.2345),{alpha,beta}); #Example kappa = 1.2345

To determine all 3 you need 3 equations.

@williamov Since you are using both 'i' and 'I' I wasn't sure what 'i' was. Changing 'i' to 'I' certainly immediately gives you a result for a.

If I understand you correctly, you don't get A(0) = 6 when n = 3, but A(0) = 33/5:

(Added: Well, I notice that you have corrected 6 to 33/5).

restart;
RA:=A(k)=1+((n-k-2)/n)*A(k+1)+(2/n)*B(k+1)+(k/n)*(A(0));
RB:= B(k)=1+((n-k-1)/n)*B(k+1)+(k/n)*A(0);
#RA holds for k = 0 .. n-2 and RB holds for k = 1 .. n-1: 
n:=3: #Example
SA:=seq(RA,k=0..n-2);
SB:=seq(RB,k=1..n-1);
solve({SA,SB},{seq(A(k),k=0..n-2),seq(B(k),k=1..n-1)});

The ode to fit the answer has y instead of y'.
eq:=diff(y(t),t,t)+y(t) = 2*t;
ics:=y(1/4*Pi) = 1/2*Pi, D(y)(1/4*Pi) = 2-sqrt(2);
sol1 := dsolve({eq, ics}, y(t), method = laplace);
expand(%);
#By "hand":
with(inttrans);
eqx:=PDEtools:-dchange({t=tau+Pi/4,y(t)=x(tau)},eq,[tau,x]);
laplace(eqx,tau,s);
eval(%,{x(0)=Pi/2,D(x)(0)=2-sqrt(2)});
isolate(%,laplace(x(tau),tau,s));
invlaplace(%,s,tau);
subs(x(tau)=y(t),tau=t-Pi/4,%);

I have revised your code.
(1) linalg is deprecated and has been replaced by LinearAlgebra
(2) I have chosen to let the B's be expressions in t and not functions. You never use them with any other argument than t.
(3) Because of the huge number q/m you will have more success using dsolve/numeric in combination with odeplot from the plots package.

restart;
with(LinearAlgebra):
with(DEtools):

SecTer := CrossProduct(<B[x], B[y], B[z]>,<diff(x(t), t), diff(y(t), t), 0>);

ProjX := diff(x(t), t, t) = -q*DotProduct(SecTer,<1, 0, 0>,conjugate=false)/m;

ProjY := diff(y(t), t, t) = -q*DotProduct(SecTer,<0, 1, 0>,conjugate=false)/m;


B[x] := 0;

B[y] := 0;

B[z] := 1 + diff(y(t),t);

q := 1.6*10^(-19); m := 9.10^(-31);

res:=dsolve({ProjX, ProjY,x(0) = 0, y(0) = 0, D(x)(0) = 0, D(y)(0) = 1},numeric,stiff=true):
plots:-odeplot(res,[x(t),y(t)],0..40);
plots:-odeplot(res,[t,y(t)],0..1e-10);
# Compare with the problems with DEplot:
DEplot([ProjX, ProjY], [x(t), y(t)], t = 0 .. 0.0004, [[x(0) = 0, y(0) = 0, (D(x))(0) = 0, (D(y))(0) = 1]],stepsize=1e-8,scene=[x,y]);

Question: Why didn't you expect solve to just express d in terms of s, a, b, and c? Should solve guess that you wanted to eliminate s?

restart;
eq1:=s=a+b+c+d;
eq2:=s=2*(a+b)*(a+c)*(b+c);
eliminate({eq1,eq2},s);
solve(%[2],d);
#Alternative:
solve(eq1,{d});
subs(eq2,%);
expand(%); #For comparison only
#But you could have eliminated e.g. c:
eliminate({eq1,eq2},c);
solve(%[2],d);

If the list is L, then

select(type,L,even);