You should not use a name and the same indexed name in the same code.

Example:
 

f := a*x+b; 
f[1] := subs(x = 1, f); 
f[2] := subs(x = 2, f);

                                                    f := a x + b
                                                  f[1] := a + b
                                                     f[2] := f

Change names in your code, for example   psi[i]  with  psi1[i]  and  Omega[i]  with  Omega1[i]

We can easily find one solution to this equation in the form of a constant function  f(x) = C :

restart;
eq:=(1/24)*exp(-8)-(1/12)*exp(-5)+(1/24)*exp(-2)-(1/24)*exp(-10)+(1/24)*exp(-9)+(1/24)*exp(-1)+diff(f(x), x, x, x, x, x)+d*(diff(f(x), x))+e*f(x)+a*(diff(f(x), x, x, x, x))+b*(diff(f(x), x, x, x))+c*(diff(f(x), x, x))-1/24;
eval(eq,f(x)=C);
f(x)=solve(%, C);
odetest(%,eq);

The initial conditions for obtaining this solution are obvious. It is also obvious that  limit(C, x=infinity)=C  for C=constant

restart;
sin(2*n*Pi) assuming posint;
sin(n*Pi) assuming integer; # More general result 

                                     0
                                     0

 minhthien2016 180 To check, you can use the procedure below, which finds the center and radius of the insphere from the coordinates of the vertices of the tetrahedron.

restart:
Center_and_Radius:=proc(P1::Vector,P2::Vector,P3::Vector,P4::Vector)
local s1, s2, s3, s4, s, V, R;
uses LinearAlgebra;
s1:=1/2*Norm(CrossProduct(P3-P2,P4-P2),2, conjugate=false):
s2:=1/2*Norm(CrossProduct(P3-P1,P4-P1),2, conjugate=false):
s3:=1/2*Norm(CrossProduct(P2-P1,P4-P1),2, conjugate=false):
s4:=1/2*Norm(CrossProduct(P2-P1,P3-P1),2, conjugate=false):
s := s1+s2+s3+s4:
V:=1/6*Determinant(<P2-P1|P3-P1|P4-P1>);
R:=3*V/(s1+s2+s3+s4);
simplify(1/s*(s1*P1+s2*P2+s3*P3+s4*P4));
convert(%, list),simplify(R);
end proc:

Your example:

Center_and_Radius(<1, 2, 3>, <-2, 8, 9>, <5, 0, 7>, <3, 4, 2>);

                                         [2, 3, 4], 1

This tetrahedron has all flat angles at the vertix  [1,2,3]  equal to 90 degrees, so you can find an equal tetrahedron in the list from my first answer.

The method above does not work in Maple 2018.2.
First, we prove the identity for integers (no problem here). If the number is not an integer, then we use the obvious fact that any not integer real number  x  can be represented in the form  x=x0+t  , x0 is an integer and  t>0, t<1 .

restart;
x:=x0+t:
ceil(x0)+floor(-x0) assuming x0::integer;
ceil(x)+floor(-x) assuming x0::integer, t>0, t<1;
is(%=0) assuming t>0,t<1;

 Earl 570 Since your equation is written in vector form, it is equivalent to a system of 3 equations with three unknown functions. Of course Maple easily solves this system. To obtain a specific solution, you need to set the initial conditions. I just formally solved these equations without delving into the physical meaning of all this:

restart;
with(LinearAlgebra):
r(t):=<r__x(t),r__y(t),r__z(t)>:
v(t):=diff(r(t),t);

Omega:=<0,0,1>: g := 9.81: theta := Pi/8:
diff(v(t), t)=~ 2/7*(Omega &x v(t)) + <0,5/7*g*sin(theta),0>;
Sys:=convert(%, list);
Sol:=dsolve(Sys);
eval(Sol,[_C1=0,_C2=0,_C3=0,_C4=0,_C5=0,_C6=0]);

We see if in the general solution we take all the undefined constants  _C  to be zero, then, in fact, motion along a straight line is obtained in the opposite direction of x-axis.

To plot something specific, you need to set the expression  g(x)  and 5 initial conditions for your 5th order equation.

Example:

restart;
g(x):=sin(x):
eq1 := diff(f(x),x$5)+2*diff(f(x),x$4)+diff(f(x),x$3)+diff(f(x),x$2)+2*diff(f(x),x)+3*f(x) = g(x);
sol:=dsolve({eq1, f(0)=0,D(f)(0)=1,(D@@2)(f)(0)=-1,(D@@3)(f)(0)=2,(D@@4)(f)(0)=-2}, numeric);
plots:-odeplot(sol, [x,f(x)], x=0..8);

If you don't like an expression with sine instead of cosine, you can easily convert sine to cosine using the  convert  command:

restart;
eqaux_2 := l__bb = L__aa0 + L__aa2 * cos(2*theta - 4/3*Pi);
eqaux_2 :=convert(eqaux_2, cos);

Here we get the expression with a cosine even simpler than the original one  ( cos(2*theta - 4/3*Pi) = cos(2*theta+2*Pi/3)  is an identity).

You can use the  plots:-display  command to place multiple graphs on one plot.

Example:

restart;
sys := {diff(x(t), t) = y(t), diff(y(t), t) = -x(t)-(1/2)*y(t)}:
A:=DEtools:-DEplot(sys, [x(t), y(t)], t = 0 .. 15, [[x(0) = 1, y(0) = 0]]):
B:=plot(-x^2, x=-0.5..1):
plots:-display(A,B);

Seems like a bug. But  assume  on a separate line works even without  simplify  and  real  (in Maple 2018.2):

restart: with(Physics):
assume(a>b[2]):
csgn(a-b[2]);  # OK

The code below finds some examples of tetrahedra with these properties. Your example is the first on the list  L :

restart:
CenterOfInSphere:=proc(P1::Vector,P2::Vector,P3::Vector,P4::Vector)
local s1, s2, s3, s4, s;
uses LinearAlgebra;
s1:=1/2*Norm(CrossProduct(P3-P2,P4-P2),2, conjugate=false):
s2:=1/2*Norm(CrossProduct(P3-P1,P4-P1),2, conjugate=false):
s3:=1/2*Norm(CrossProduct(P2-P1,P4-P1),2, conjugate=false):
s4:=1/2*Norm(CrossProduct(P2-P1,P3-P1),2, conjugate=false):
s := s1+s2+s3+s4:
simplify(1/s*(s1*P1+s2*P2+s3*P3+s4*P4));
convert(%, list);
end proc:

N:=20:
k:=0:
for a from 1 to N do
for b from a to N do
for c from b to N do
C:=CenterOfInSphere(<0, 0, 0>, <a, 0, 0>, <0, b, 0>, <0, 0, c>);
R:=C[1];
if R::integer then k:=k+1; L[k]:=[[a,b,c],R] fi;
od: od: od:
L:=convert(L, list);

L := [[[3, 5, 14], 1], [[3, 6, 9], 1], [[4, 4, 8], 1], [[6, 12, 18], 2], [[6, 13, 16], 2], [[8, 8, 16], 2], [[8, 10, 11], 2], [[12, 14, 18], 3]]

You can increase the  N  number to get more solutions.

Use  simplify  with side relation for this:

restart;
simplify(exp(-a), {exp(a)=p});

Also

f:=exp(-a)*(b+c)*1/k;
simplify(f, {exp(a)*k=p});

The calculations on my machine (Maple 2018) took only 23 seconds. Found 64 solutions:

restart;
ts:=time():
k:=0:
for x from 3 to 30 do
for y from 2 to x-1 do
for z from 1 to 30 do
for a from 3 to 10 do
for b from 2 to a-1 do
for n from 2 to 10 do
for m from 2 to n do
if igcd(m,n)=1 and igcd(a,b,m)=1 then 
if 1/x+1/y=1/z and a/x+b/y=m/n then k:=k+1; L[k]:=[x,y,z,a,b,m,n]
fi; fi;
od: od: od: od: od: od: od:
L:=convert(L, list);
nops(L);
time()-ts;

Equality  x+y=1  is possible only if  y=1-x . But then the contradiction disappears:

diff(x+(1-x)=1, x);

                          0 = 0

For other values of  y , there is no any equality. Therefore, it is not surprising that the derivatives do not coincide.

General conclusion: if some equality  f(x)=g(x)  is an identity (that is, it is true for all values of the variable  x) then  f(x)=g(x) => diff(f(x),x)= diff(g(x),x) . If equality is not an identity, then both cases are possible  diff(f(x),x)= diff(g(x),x)  as well as diff(f(x),x)<>diff(g(x),x) .

series(1/(1-x^2), x=1);

Or

numapprox:-laurent(1/(1-x^2), x = 1)

This expansion is hold in the ring  0<|x-1|<2  containing no singular points  x=1  and  x=-1 .