We prove a comprehensive solution theory using tools from functional analysis, show corresponding variational formulations, and present functional a posteriori error estimates for general linear first order systems of type A₂x-f,A^* ₁x=g, for two densely defined and closed (possibly unbounded) linear operators A₁ and A₂ having the complex property A₂A₁=0. As a prototypical application we will discuss the system of electro-magneto statics in 3D with mixed tangential and normal boundary conditions rotE-F - div єE=g, Our theory covers a lot more applications in 2D, 3D, and ND, such as general differential forms and all kind of systems arising, e.g., in general relativity, biharmonic problems, Stokes equations, or linear elasticity, to mention just a few, for example dE=F, Rot_SM=F, Div_TT=F, Rot Rot^T _S S=F, -ᵟєE=G, div Div_sєM=G, sym Rot_TєT=G,-Div_S єS=G,all with possibly mixed boundary conditions of generalized tangential and normal type. Second order systems of types A^* ₂A₂x=f,A^* ₂A₂x=f, A^* ₁x=g, A₁A^* ₁x=g, will be considered as well using the same techniques.