We need to see the set of ODEs.

@geri23 From your description it seems to me that the system would be something like:

C1*diff(T1(t),t)=U12*(T2(t) -T1(t))+ U01*(Tout -T1(t))+ H1(t);
C2*diff(T2(t),t)=U12*(T1(t) -T2(t))+ U02*(Tout -T2(t))+H2(t);


where Tout is the outside temperature and the constants Uij > 0, Ci > 0. You didn't say what H1(t) or H2(t) were, but they must be heating sources inside the two rooms. In order to solve the system, you need information about those:  Are they in fact known constants?
To solve the system you need the initial values of T1 and T2.
dsolve should be able to solve this linear system with no problem (certainly no problem if H1 and H2 are constants).

@Kivo69 The integral in this problem won't go away by differentiation.
Maybe replacing the integrand with a series expansion about xi would at least initially be interesting:

IntegrationTools:-GetIntegrand(eq3);
series(%,zeta=xi,2);
convert(%,polynom);
Int(%,zeta=0..1);
value(%);
eq3a:=evalindets(eq3,specfunc(anything,Int),x->u);
#I haven't had time to look further

Are i and j indices or just names?
If they are indices, say i = 1 .. n , j = 1 .. n, then you could clarify the question by writing out all the equations when n = 2. Can i be equal to j?
If i and j are just names and you only have one equation then you need to know what Tj(t) is (and Hi,Uij,Ci of course).

@KatePirs In one of your by now numerous postings of this question it has been pointed out (I think by Carl Love) that your x-range might as well be considerably smaller. Anyway, try this:

IBC2 := {ui(-10, t) = 0, ui(10, t) = 0, ui(x, 0) = 0, ur(-10, t) = 0, ur(10, t) = 0, ur(x, 0) = 1/cosh(x)};
sol := pdsolve({l1, l2}, IBC2, funcs, numeric, time = t);
sol:-plot([ui,ur],t=10);
sol:-animate([ui,ur],t=500,frames=100);
sol:-plot(ur,t=0,x=-10..10,numpoints=500);
plot(1/cosh(x),x=-10..10);


The similarity is striking.

@KatePirs In one of your by now numerous postings of this question it has been pointed out (I think by Carl Love) that your x-range might as well be considerably smaller. Anyway, try this:

IBC2 := {ui(-10, t) = 0, ui(10, t) = 0, ui(x, 0) = 0, ur(-10, t) = 0, ur(10, t) = 0, ur(x, 0) = 1/cosh(x)};
sol := pdsolve({l1, l2}, IBC2, funcs, numeric, time = t);
sol:-plot([ui,ur],t=10);
sol:-animate([ui,ur],t=500,frames=100);
sol:-plot(ur,t=0,x=-10..10,numpoints=500);
plot(1/cosh(x),x=-10..10);


The similarity is striking.

@helix  Something like this:

with(plots):
p1:=animate(arrow,[<cos(tend),sin(tend),0.9*tend>,<-cos(tend),-sin(tend),0>,color=red],tend=0..50,axes=boxed,frames=150);
alpha:=1:
p2:=animate(textplot3d,[[[-1.8,0,alpha*tend,piecewise(tend<=12,"",tend<=25,"position","")]]],tend=0..50,axes=boxed,frames=150);
display(p1,p2);

@helix  Something like this:

with(plots):
p1:=animate(arrow,[<cos(tend),sin(tend),0.9*tend>,<-cos(tend),-sin(tend),0>,color=red],tend=0..50,axes=boxed,frames=150);
alpha:=1:
p2:=animate(textplot3d,[[[-1.8,0,alpha*tend,piecewise(tend<=12,"",tend<=25,"position","")]]],tend=0..50,axes=boxed,frames=150);
display(p1,p2);

There are various ways. But we shall need the details, i.e. all the lines of code.

@KatePirs
I don't get that error message, but after supplying values for a and m I get complaints about initial/boundary conditions.
You have boundary conditions, but no initial conditions!

The help page for pdsolve,numeric says:

"The pdsolve(PDEsys,conditions,numeric,vars,options) command returns a module that can be used to compute numerical solutions for time-based PDE systems over a fixed finite 1-space interval."


That must mean that initial conditions must be given also.

@KatePirs
I don't get that error message, but after supplying values for a and m I get complaints about initial/boundary conditions.
You have boundary conditions, but no initial conditions!

The help page for pdsolve,numeric says:

"The pdsolve(PDEsys,conditions,numeric,vars,options) command returns a module that can be used to compute numerical solutions for time-based PDE systems over a fixed finite 1-space interval."


That must mean that initial conditions must be given also.

@Axel Vogt Actually the worksheet was saved with the data. If you execute the assignment to stiff you will see that it is a symmetric tridiagonal matrix with diagonal elements = 2 and the other -1, as would have been produced by:
LinearAlgebra:-BandMatrix([-1,2,-1], 1, 21);

It would be more interesting and useful if you gave us the 21x21 matrix. Maybe upload a worksheet.

Works for me too in 16.02.

Works for me too in 16.02.